Unravelling raked linear dunes to explain the coexistence of bedforms in complex dunefields

Raked linear dunes keep a constant orientation for considerable distances with a marked asymmetry between a periodic pattern of semi-crescentic structures on one side and a continuous slope on the other. Here we show that this shape is associated with a steady-state dune type arising from the coexistence of two dune growth mechanisms. Primary ridges elongate in the direction of the resultant sand flux. Semi-crescentic structures result from the development of superimposed dunes growing perpendicularly to the maximum gross bedform-normal transport. In the particular case of raked linear dunes, these two mechanisms produces primary and secondary ridges with similar height but with different orientations, which are oblique to each other. The raked pattern develops preferentially on the leeward side of the primary ridges according to the direction of propagation of the superimposed bedforms. As shown by numerical modelling, raked linear dunes occur where both these oblique orientations and dynamics are met.

R aked linear dunes display asymmetric morphologies with long primary ridges breaking down only on one side into periodic secondary ridges with a perpendicular alignment to the main crest and a semi-crescentic shape. This unique pattern provides an intriguing case study to begin to examine how primary and secondary bedforms with similar heights but different trends can coexist in the same wind regime.
A major question in arid zone geomorphology is to determine to what level active dunefields are in dynamic equilibrium with the current wind regime or the result of longer-term processes like climate change 1 . This is a challenging task given the broad spectrum of observed dune shapes, which are primarily controlled by sand availability and wind directional variability 2,3 . If crescentic barchans (low sand availability, unidirectional wind) and star dunes (high sand availability, multidirectional winds) occur under two extreme and opposite sets of conditions, linear dunes are often sub-divided into different categories according to their shape 4,5 (simple, compound and complex) or orientation with respect to the resultant sand flux direction 6,7 (transverse, oblique and longitudinal). Unfortunately, all these classifications fail to take account of dynamics and, without independent chronological data [8][9][10] , it is difficult to interpret or reconstruct the dunefield history.
Rather than a single growth mechanism, it was recently shown 3,11 that dunes can either extend in the direction of the resultant sand flux in zones of low-sediment supply 12 (that is, the fingering mode) or grow perpendicularly to the maximum gross bedform-normal transport 13 in conditions of high-sediment supply (that is, the bed instability mode). These two distinct dune growth mechanisms determine both dune shape and orientation. Taking into account the feedback of dune aspectratio on the magnitude of the flow 14 , these orientations can be analytically derived from the wind regime to be compared with observations in the field. Most importantly for our present purpose, the coexistence of the two dune growth mechanisms can theoretically lead to the development of different types of steadystate bedforms through spatial changes in sand availability or the emergence of secondary structures 15 . Given the apparent complexity of most dunefields on Earth, none of these interacting steady-state dune patterns have been identified and quantified so far in modern sand seas. It is an important issue because the coexistence of different bedform alignments has also been associated with changes in wind regime 16 .
Raked linear dunes have been recognized for the first time in the Kumtagh desert 17 (Xinjiang province, China), where they cover an area of more than 2.4 Â 10 4 km 2 . This remote and arid region of China has been the subject of intensive investigation in recent years [18][19][20][21][22] . These studies provide unique sets of data, but none of them have identified the mechanism for dune formation and the subsequent dune morphodynamics. Instead, most of them have concentrated on zibars and other grain-size segregation patterns 23 , which are obvious in interdune areas.
Considering the traditional classification, raked linear dunes seem to be a specific case of compound/complex linear dunes. Their most distinguishable feature is a lateral asymmetry in the distribution of secondary ridges. At the dunefield scale, primary ridges display all the characteristics of linear dunes. They are relatively straight, extending over distances of kilometres parallel to one another with a low ratio of dune to interdune areas 24 . On the top of them, secondary ridges almost perpendicular to the main dune alignment develop only on one side of the dune. These secondary ridges have approximately the same height as the primary ridges and they resemble half crescent, so that the entire pattern is similar to a train of barchans connected to each other along the same arm.
Here we study the morphodynamics of raked linear dunes in the Kumtagh desert from two satellite images taken at a time interval of 8 years. Dune shapes, orientations and dynamics (elongation and migration) are compared with theoretical predictions derived from local wind data to show that raked linear dunes are a steady-state dune type in dynamic equilibrium with the modern wind regime. Then, we numerically and analytically investigate the conditions for the emergence of this dune type to demonstrate that the raked pattern result from the coexistence of two dune growth mechanisms.

Results
Morphodynamics of raked linear dunes. Figure 1a,b shows aerial views of raked linear dunes in the Kumtagh desert (see also Supplementary Figs 1-3). At different sites in interdune areas of the Kumtagh desert, two wind towers of 2 m height have recorded the mean wind speed and direction every 15 min during 2 years from 2007 to 2009. Figure 1c shows the wind roses and the distributions of sand flux direction. This eastern end of the Tarim basin is exposed to a trimodal wind regime with a widely spread dominant wind direction from the north and two secondary peaks associated with easterly and westerly winds. Dune orientations predicted from these wind data are a north-south alignment for bedforms growing in the bed instability mode and an extension in the southwesterly direction for dunes growing in the fingering mode ( Fig. 1d and Table 1 in 'Methods' section). These orientations are in general agreement with the alignments of the primary and secondary ridges observed in satellite images (Fig. 1e). These findings suggest that raked linear dunes can form by elongation and that secondary ridges can result from the development of superimposed bedforms growing perpendicularly to the maximum gross bedform-normal transport. In the framework of the two dune growth mechanisms, the two dune orientations form an oblique angle Da ¼ 70 ± 3°, which leads to an angle Da Q ¼ 5 ± 2°b etween the resultant sand flux at their crests because of the speed-up effect in multidirectional wind regimes 15 (see 'Methods' section). Then, as the superimposed bedforms nucleate and grow on the primary ridges, they migrate from one side of the linear dunes to the other. According to this migration, the regular pattern of secondary ridges is observed on the leeward slope of the primary ridges (that is, the north-western side of the linear dunes in Fig. 1). Hence, the complex asymmetric pattern of raked linear dunes may naturally arise from the coexistence of two dune growth mechanisms with both oblique orientations and dynamics.
Independently of the winds, such a scenario can be quantitatively investigated from dune morphodynamics using field data and satellite images from the Kumtagh desert ( Fig. 2a,b). Figure 2c shows that the wavelength of the semicrescentic secondary bedforms is linearly related to the height of the linear dunes. This relationship is similar to that observed for longitudinal dunes breaking in barchans under unidirectional wind regime [25][26][27] , suggesting that the same mechanism is responsible for the singular asymmetry of raked linear dunes. In addition, the widths W B of the semi-crescentic structures exhibit a normal distribution with a mean value around 37 m and a s.d. close to 3 m across the entire dunefield (Fig. 2d). The linear dunes reach a homogeneous width W F with a mean value of 26 m (Fig. 2e) away from the growing tips, which may change shape over time (see 'Methods' section). Using a pair of Google Earth images acquired in 2008 and 2014, we estimate the elongation rate e ¼ 13.2 ± 3 m per year at the tip of linear dunes and the migration rate c ¼ 6±1 m per year of the secondary ridges along the same direction (Fig. 2f). Despite more variation in the elongation rate distribution, which can be explained by the loss of sediment at the dune tips during periods of constant wind orientation, Fig. 2g shows that linear dunes elongate more than two times faster than the secondary ridges migrate.
Sand fluxes on raked linear dunes. Using the morphometric parameters of raked linear dunes, the resultant sand fluxes Q F and Q B at the crests of primary and secondary ridges can be directly derived from the equation of conservation of mass over the entire cycle of wind reorientation. Considering that the dunes elongate or propagate without changing shape, the relations are e ¼ Q F /H F and c/cos(Da Q ) ¼ Q B /H B , where cos(Da Q ) is a correction to account for the migration direction of the superimposed bedforms in the bed instability mode and where H F and H B are the heights at the crests of primary and secondary bedforms, respectively. These heights H F and H B are derived from the measurements of the width W F of linear dunes (Fig. 2b,f) and from the width W B of half-crescentic structures (Fig. 2b,e) using the rela- . Here we take y F ¼ 15°a nd y B ¼ 11°, which are the usual reported values for the transverse aspect-ratio of linear 12 and barchan 28 dunes. Figure 3 shows the comparison between the sand fluxes computed from wind data and dune morphodynamics. The close agreement between these two independent estimates indicate that a hierarchy of bedforms can be used to provide independent quantitative constraints on sediment transport.
Conditions for steady-state raked linear dunes. Together, sand flux estimates (intensity and direction) and comparisons between predicted and observed dune orientations suggest that the raked linear dunes in the Kumtagh desert have reached a steady state,  (d) Predicted orientations for dunes growing by extension (fingering mode, black arrows) and perpendicularly to the maximum gross bedform-normal transport (bed instability mode, red arrows). Green arrows show the predicted resultant sand flux at the crest of dunes growing in the bed instability mode. By definition, this is also the propagation direction of these dunes. (e) Comparison between predicted and observed dune orientation at two different times.
which is fully consistent with the recorded wind regime. Nevertheless, such a regular pattern has never been reported before in numerical and laboratory experiments investigating dune shape 11,25,[29][30][31][32] . The main reason is that all these studies have only concentrated on unidirectional and bidirectional wind regimes that seem to never satisfy the conditions for the development and stability of raked linear dunes (see 'Methods' section). From our field-based results and systematic exploration of dune shape in bidirectional wind regimes 11 , we infer that raked linear dunes can only be observed in zones of low sand availability if three conditions are met (see 'Methods' section and equations (2)(3)(4)(5)(6)(7)(8)(9)(10) for the derivation of all variables from wind data). First, the ratio s F /s I between the dune-height growth rates of the fingering and the bed instability modes has to be large enough to promote elongation. Second, the angle Da between the two dune orientations should be oblique, wide enough for the two modes to be distinguished and yet narrow enough to avoid segmentation of the primary ridges (especially for decreasing s F /s I -value). Finally, the angle Da Q between the migration direction of the superimposed bedforms and the elongation  direction should be wide enough to preserve the continuity of primary ridges and yet narrow enough for the secondary ridges and the asymmetric pattern to fully develop.
To test these hypotheses, we use a cellular automaton dune model that accounts for feedback mechanisms between the flow and the evolving bed topography 11,12,15,[33][34][35] . Using a tridirectional wind regime similar to that observed in the field, simulations are run for s F /s I ¼ 0.75, Da ¼ 42°and Da Q ¼ 9.3°. After a few cycles of wind reorientation a linear dune rapidly extends on the non-erodible bed away from the fixed source of sediment (Fig. 4). As it continues to elongate, superimposed bedforms develop and grow preferentially on the leeward side of the linear dunes according to their direction of propagation. Then, a systematic asymmetry starts to develop despite the constant loss of sediment at the tips of the linear dune and semicrescentic structures. Over longer times, the raked pattern extends across the entire domain, keeping a constant shape characterized by regularly spaced secondary ridges of constant height and width propagating only on one side of a well-established linear dune. The tip region left appart, the superimposed dunes do not escape the finger dunes, but migrate along it. It is worth noting that the raked linear dunefield in the Kumtagh desert does not exhibit isolated barchan either. As in the numerical simulations, this suggest that the mass balance between the primary and the secondary ridges has reached an equilibrium.

Discussion
Numerous observations of differing alignments between main and superimposed bedforms have previously been observed in both modern aeolian and subaqueous bedforms as well as in their ancient deposits. For example, Rubin and Hunter 36 and Rubin 37,38 present examples of deposit records where the two sets of bedforms coexist under (inferred) steady-state conditions. These previous studies attributed differences in orientations not to changes in flow direction through time but rather to differences between the flow generating the main bedforms (external to the dune field) and the flow producing the superimposed dunes (flow within the internal boundary layer).
Here we find that, in multidirectional wind regimes, there is no need for secondary flow or deflection of lee-side flow to give rise to the coexistence of bedforms with similar heights and different trends. There is no such ingredient in our numerical model. Instead, considering specific boundary conditions and sediment properties 23 , all the independent variables required to produce the different types of dunes are incorporated into the function of flow directionality and intensity. This approach is especially justified when it comes to study primary and secondary ridges with the same height.
Our results show that raked linear dunes can be considered as a new class of composite bedforms, just as star dunes 15 . According to the usual classifications, raked linear dunes can be considered as linear dunes asymmetrically indented by a train of superimposed oblique dunes. Such a hierarchy of bedforms is possible because primary and secondary ridges are associated with different dune growth mechanisms with specific orientation and dynamics. The main dunes are in supply-limited conditions and, as a result they elongate in the direction of the resultant sand flux (fingering mode). The superimposed dunes arise from the instability of a sand bed 39 because the supply of sand from the main longitudinal dunes is essentially unlimited. Nevertheless, the simultaneous manifestation of these two modes of dune orientation occur only in limited angular ranges between crest alignments and their propagation or extension directions. As shown by numerical and analytical modelling, a tridirectional wind regime similar to that observed in the field is needed to fulfil these conditions. Unidirectional or bidirectional wind regimes cannot 11 (see 'Methods' section), which is consistent with the rareness of this dune type.
From the ratio between the elongation rate (E13 m per year) and the length of the raked linear dunefield (E100 km), the minimum time for the formation of the Kumtagh desert dunes under the present wind conditions is at least 5,500 years. Such a time scale and the principal wind directions are in good agreement with the presence and the alignment of yardangs at the northeastern end of the dunefield ( Supplementary Fig. 4), indicating that the current wind regime may have prevailed since the Holocene thermal maximum.
More generally, this study illustrates how the combination of two dune growth mechanisms can generate new types of steadystate bedforms in dynamic equilibrium with a given multidirectional wind regime. Complex dune shapes with multiple crest orientations may develop from the coexistence of bedforms with the same 15,40 or different 3 growth mechanisms. Considering the possible number of combinations, it is impossible at this stage to generalize the dynamical behaviour observed in the Kumtagh to other environments where raked patterns have been observed. For example, there are similar dune shapes in Mali (Fig. 5a), the Namib Sand Sea 41 (Fig. 5b) and Saudi Arabia 42 (Fig. 5c). Nevertheless, these dunes are one order of magnitude larger than those in the Kumtagh desert and they develop in the middle of sand flow paths from other dune types over long time scales (4100 years) for which there is no wind data. As a consequence, their origin and stability need to be investigated further before any conclusion can be drawn. Indeed, they may also reflect transitions in dune type across space or over time 25,26 . Another interesting systematic feature to be explored is the raked dune patterns on Titan 43 , the largest moon of Saturn, where trimodal wind regimes have been proposed 44,45 to explain the eastward elongation of linear dune at a global scale.
Interacting bedforms have been a challenge for decades to address major issues on the dynamics of dunefields, their origin and evolution in presence of various environmental variables [46][47][48][49] . The methodology based on the two dune growth mechanisms applies generally to dune landscapes in which different trends coexist. Here we show how it can be used to

Methods
Saturated flux on a flat sand bed. Using the wind data collected in the Kumtagh desert (Fig. 1c), Table 1 shows the predicted sand flux on a flat sand bed and a set of variables relevant for dune morphodynamics: orientation, sand flux at the crest, migration direction and dune-height growth rate. Continuous wind measurements in the field provide the wind speed u i and directionx i at different times t i , iA[1; N].
For each time step i, we calculate the shear velocity where z ¼ 2 m is the height at which the wind velocity u i has been measured, z 0 ¼ 10 À 3 m the characteristic surface roughness and k ¼ 0.4 the von-Kármán constant. The threshold shear velocity value for motion inception can be determined using the formula calibrated by Iversen and Rasmussen 51 Using the gravitational acceleration g ¼ 9.81 m s À 2 , the grain to fluid density ratio r s =r f ' 2:05Â10 3 and the grain diameter d ¼ 180 mm, we find u c ¼ 0.19 m s À 1 , which corresponds to a threshold wind speed of 3.6 m s À 1 two metres above the ground. For each time step i, the saturated sand flux Q i ! on a flat sand bed can be calculated from the relationship proposed by Ungar and Haff 52 and calibrated by Durán 53 Here the sand flux takes into account a dune compactness of 0.6. From the individual saturated sand flux vectors Q i ! , we estimate the norm of the mean sand flux on a flat erodible bed, also called the resultant drift potential: (e) A snapshot of the steady-state raked linear dune. Insets show the sand flux roses and the predicted orientations of dunes growing by extension (fingering mode, black arrow) and perpendicularly to the maximum gross bedform-normal transport (bed instability mode, red arrows). The green arrow shows the resultant sand flux at the crest of dune growing in the bed instability mode. By definition, this is also the migration direction of these dunes. The asymmetric dune pattern results from the obliquity between the two dune orientations and the oblique propagation of secondary bedforms.
where dt i ¼t i À t i À 1 : This quantity is strongly dependent on the function of wind directionality. We also calculate the drift potential, which is the average of individual norms and represent the mean volume of sand per unit width and time that have been transported through a vertical plane: Averaged over the entire time period, this quantity does not take into account the orientation of the sand fluxes 54 . The ratio RDP/DP is a non-dimensional parameter, which is often used to characterize the directional variability of the wind regimes 55,56 : RDP/DP-1 indicates that sediment transport tends to be unidirectional; RDP/DP-0 indicates that most of the transport components cancel each other. Finally, RDD is the resultant drift direction, that is, the direction of Sand flux at the crest of dunes. A positive topography accelerates the wind, so that the sand flux over a dune depends on the dune shape. For 2D turbulent flows over low hills, Jackson and Hunt 57 show analytically that the increase of wind velocity at the top of the hill, the so-called speed-up factor, is approximately proportional to the hump aspect ratio. Hence, at the first order of the dune aspect ratio and neglecting the transport threshold, the sand flux Q c i ! at the crest of the dune takes the following form where a is the orientation angle of the linear dune, Q i ! the sand flux on a flat sand bed, y i the orientation angle of this flux vector and g the flux-up ratio: where W is the width of the dune, H its height and b a dimensionless coefficient that accounts for all the other physical ingredients (for example, roughness) that affect the speed-up.
Dune orientation a I in the bed instability mode. In the bed instability mode, dunes grow in height from the erosion of the sand bed in the interdune areas. All winds contribute to the growth of the dune height. Thus, linear bedforms develop perpendicularly to the maximum gross bedform-normal transport 13 .
Considering the orientation angles y i of fluxes Q c i ! , we calculate Q > (a), the total sand flux perpendicular to the crest for all possible crest orientations aA[0; p]. Then, we identify the maximum value of Q > (a) that corresponds to the most probable crest orientation a I of dunes in the bed instability mode. Note that this procedure is the same as in Rubin and Hunter 13 (that is, the gross bedform-normal transport rule), except that we take into account the increase of the sand fluxes at the crest of dunes (that is, ga0). As detailed in Courrech du Pont et al. 3 and Gao et al. 11 , it may significantly change the predictions of dune orientations.
Dune orientation a F in the fingering mode. In the fingering mode, dunes extend in the direction of the mean sand flux. Hence, the orientation a F of finger dunes is the one for which the direction of the mean sand flux at the crest (that is, the time average of equation (6)) aligns with dune orientation.
In practice, we calculate Q > (a) and Q || (a), the total sand flux perpendicular and parallel to the crest for all possible crest orientations aA[0; 2p]. Then, we select the orientation a F for which the sediment flux perpendicular to the crest vanishes (that is, Q > (a) ¼ 0) and for which the flux parallel to the dune is positive (that is, Q || (a)40). If more than one solution exists, as it is the case for star dunes 15 , we look for the angle at which the Q || -value is maximum. By definition, when there is no feedback of topography on the flow (that is, g ¼ 0 in equation (6)), the orientation of the linear fingering mode a F is given by the resultant sand transport direction (also called the RDD). When the wind speed-up is taken into account, the dune orientation depends on the g-value; g ¼ 1.6 gives reasonable estimates of dune orientation 3,11 .
From the predicted crest orientations a I and a F in the bed instability and the fingering modes, we calculate the angle between the two bedform alignments (DaA[0; p/2]).
Sand flux at the crest of dunes in the bed instability and fingering modes.
Because the apparent dune aspect-ratio (that is, the aspect ratio seen by the wind) determines the increase in wind speed at the top of the dune, the magnitude and (c) Namib desert, Namibia (24°15 0 S, 14°59 0 E). Despite similar raked patterns, we cannot conclude at this stage that these dunes are in a steady state or in equilibrium with the current wind regime. Given their size and location along major sand flow paths, they can also result from transition in dune shape across space and over time.
the orientation of the time averaged sand flux at the dune crest Q I;F f g ! is a function of the dune orientation a {I,F} : where y i is the orientation angle of the flux Q i ! on a flat sand bed. Then, Da Q is the angle between the direction of the resultant sand flux at the crest of dunes in the bed instability and the fingering modes (Da Q A[0; p]).
Dune-height growth rate. All sand fluxes perpendicular to the crest can contribute to dune growth. Considering the dune orientation a {I,F} , we calculate the characteristic (growth) rate s {I,F} to build up a linear dune of height H and width W in the bed instability or the fingering mode 3,11 : Characterization of the wind regime. We use an Expectation-Maximization algorithm to fit the flux orientation distribution by a Gaussian mixture model. Thus, we replace the real flux data by a limited number n Y of normal distributions characterized by a mean orientation Y i , a standard variation s i and a weight w i with i ¼ {1, y, n Y }. Considering only time periods during which the wind velocity is above a critical value for sediment transport (that is, u * 4u c , see equation (3)), we assume that the probability distribution function of sand flux orientation Y may be described by a sum of normal distributions: The Expectation-Maximization algorithm is a natural generalization of maximum likelihood estimation to the incomplete data case. Basically, this is an iterative scheme that includes two different steps. Starting from initial guesses for the parameters w i , s i and Y i , the (first) Expectation-step is to compute a probability distribution over possible completions. In the (second) Maximization-step, new parameters are determined using the current completions. These steps are repeated until convergence. Elongation rate and mean width of linear dunes. To estimate the elongation and the mean width of linear dunes from Google Earth images, we isolate the contours of the growing tip at two different times ( Supplementary Fig. 5a). From the surface area covered by individual dunes, we compute dune orientation and locate the dune tip at each time. Then, we measure dune width perpendicularly to this orientation at different positions along the linear dune. Supplementary Fig. 5b-e shows the variation of dune width with respect to the distance from the latest dune tip. We measure the elongation using the distance between the dune tips at two different times. We estimate the mean dune width as the homogeneous width reached along the dune body away from the growing tip. From these measurements, we observe that the shape of the growing tip changes. This is not surprising given the function of wind directionality. For this reason, we could have measured elongation using the positions at which dunes reach their homogeneous width instead of the dune tips. Nevertheless, this procedure gives similar results but more uncertainty. One should also note that the sand fluxes derived from these elongation rates are likely to be underestimated because of sand loss at the dune tips.
No raked linear dunes under bidirectional wind regime. Raked linear dunes have never been reported before in laboratory or numerical experiments investigating dune shape 25,[29][30][31][32] , even when the two dune growth mechanisms have been under consideration 3,11 . All these experiments have concentrated only on unidirectional or bidirectional wind regimes. There is clearly three dominant wind directions in the Kumtagh desert. Such a trimodal wind regime generates specific conditions, which are never met together in a bidirectional wind regime.
In the parameter space {y, N} of bidirectional wind regimes (y and N are the angle and the transport ratio between the two winds, respectively), {Da, s F /s I , Da Q }-values can be computed analytically 11 (Supplementary Fig. 6a-c). For each variable, we report the range of values derived from the wind data in the Kumtagh desert to the parameter space {y, N} of bidirectional wind regimes. Thus, we delineate different regions of the parameter space {y, N}. When we compare these regions, we observe that there is no overlap indicating that the specific conditions predicted from the recorded wind data cannot be reproduced by bidirectional wind regimes ( Supplementary Fig. 6d-f). Therefore, there is a general incompatibility between the three variables and, even when they are considered two by two, regions of overlaps are small in size. Hence, bidirectional wind regimes cannot provide the {Da, s F /s I , Da Q }-values met in the raked linear dunefield of the Kumtagh desert.
Data availability. The data that support the findings of this study are available from the corresponding author upon reasonable request. The numerical dune model has been built using the Real-Space Cellular Automaton Laboratory 35 (ReSCAL), a free software under the GNU general public licence. The source codes can be downloaded from http://www.ipgp.fr/rescal.