title: Spatial and temporal development of incipient dunes
authors: Cyril Gadal Clément Narteau Ryan C. Ewing Andrew Gunn Douglas Jerolmack Bruno Andreotti Philippe Claudin
text: 
  
  
  
  
  
  
  
  
  C. Gadal1, C. Narteau1, R.C. Ewing2, A. Gunn3, D. Jerolmack3,4, B. Andreotti5, P. Claudin6
  
  1Institut de physique du globe de Paris, Université de Paris, CNRS, Paris, France
  2Department of Geology and Geophysics, Texas A&M University, College Station, USA
  3Department of Earth and Environmental Science, University of Pennsylvania, Philadelphia, USA
  4Department of Mechanical Engineering and Applied Mechanics, University of Pennsylvania, Philadelphia, USA
  5Laboratoire de Physique, ENS - PSL Research University, Université de Paris, CNRS, Sorbonne Université, Paris, France
  6Physique et Mécanique des Milieux Hétérogènes, ESPCI Paris - PSL Research University, Université de Paris, CNRS, Sorbonne Université, Paris, France
  
  
  Cyril Gadalgadal@ipgp.fr
  
  
  An edited version of this paper was published by AGU. Copyright 2020 American
  Geophysical Union: Gadal, C., Narteau, C., Ewing, R. C., Gunn, A., Jerolmack, D., Andreotti, B., & Claudin, P. (2020). Spatial and temporal development of incipient dunes. Geophysical Research Letters, 47, e2020GL088919. https://doi.org/10.1029/2020GL088919.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
    * Length and time scales of dune formation set the pace of aeolian landscape dynamics
  
    * We provide significant statistics on the temporal and spatial growth of incipient sand dunes extracted from extensive field observations
  
    * Data provide a field validation of dune instability theory, and introduce the distance of dune growth as an important length scale
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  In zones of loose sand, wind-blown sand dunes emerge due the linear instability of a flat
  sedimentary bed. This instability has been studied in experiments and numerical models but rarely in the
  field, due to the large time and length scales involved. We examine dune formation at the upwind margin
  of the White Sands Dune Field in New Mexico (USA), using 4 years of lidar topographic data to follow
  the spatial and temporal development of incipient dunes. Data quantify dune wavelength, growth rate, and
  propagation velocity and also the characteristic length scale associated with the growth process. We show
  that all these measurements are in quantitative agreement with predictions from linear stability analysis.
  This validation makes it possible to use the theory to reliably interpret dune-pattern characteristics and
  provide quantitative constraints on associated wind regimes and sediment properties, where direct local
  measurements are not available or feasible.
  
  
  
  
  § PLAIN LANGUAGE SUMMARY
  
  
  
  
  Dunes are the solar system's ubiquitous landform, arising wherever
  wind blows over a loose sand bed. An aerodynamic theory for dune formation, which connects grain-scale
  movement to emergent dune pattern, has been developed for idealized scenarios. Yet this model has
  never been directly tested in nature, because of the complexities in observing dune formation at the initial
  stage. Here we report extensive topographic observations of the initiation, growth, and migration of
  real-world sand dunes. Moreover, we find a surprisingly precise agreement with the idealized aerodynamic
  theory. This robust confirmation of the theory for dune formation means that we may estimate wind
  conditions in remote places, including other planets, with confidence.
  
  
  
  
  
  
  
  § INTRODUCTION
  
  
  
  The development of sand dunes, from incipient to mature bedforms, and their evolution, through interaction and coarsening processes, involve characteristic time and length scales that relate to elementary mechanisms of hydrodynamics and sediment transport <cit.>. Over loose granular beds, bedform emergence is driven by a hydrodynamic instability induced by the interaction between the sand bed, flow, and sediment transport <cit.>. On the upstream side of a bump, erosion takes place as the flow accelerates. Reciprocally, the flow slows down on the downstream side where deposition occurs. However, the transition between erosion and deposition zones, associated with the location of the maximum of the sediment flux, does not necessarily take place at the crest of the bump. Spontaneous growth of such a bump–that is, instability–can therefore occur if its crest is located in the deposition zone <cit.>. The streamwise offset between topography and sediment flux has two contributions <cit.>. First, a hydrodynamic destabilization originates from the coupling between flow inertia and dissipation, which results in a maxi-
  mum basal fluid shear stress located upstream of the crest <cit.>. Second, the sand flux needs a characteristic length, called the saturation length, to adapt to a spatial change in shear stress <cit.>. This results in a stabilizing downstream lag of the maximum sand flux with respect to the maximum of the
  shear stress. These balancing processes give rise to the development and propagation of sand waves at a specific wavelength and propagation speed, associated with the most unstable mode of the pattern, with crests perpendicular to the dominant wind direction.
  
  The early stage of growth and development of sedimentary ripples and dunes has been theoretically studied with linear stability analyses of coupled transport and hydrodynamic equations <cit.>. These analyses predict the incipient pattern wavelength, propagation velocity, and growth rate as functions of model parameters, which encode the wind and grain characteristics. For the aeolian case in particular, the dune wavelength has been shown to be proportional to the saturation length. However, measuring the bed elevation together with sediment and fluid transport is difficult, thereby making the direct comparison between theory and field or experimental data rather challenging.
  
  The aerodynamic and sediment transport responses have been investigated independently of each other, and separate measurements of the saturation length and the upwind shift of the shear stress have been carried out, in the field and in wind tunnel experiments <cit.>. In contrast, few field studies addressing the early stage of aeolian dune growth are available in the literature <cit.>. First, in situ monitoring of the evolution of small amplitude bedforms is difficult due to the involved length and time scales (tens to hundreds of meters, days to months). Second, inherent wind variability–even in overall unidirectional dune fields–makes application of the theory challenging. Emergence of subaqueous sand ripples has also been experimentally investigated <cit.>, and more generally the quantification of sedimentary bedforms in different environmental–including extraterrestrial–conditions in relation to hydrodynamics and sediment transport remains an active current subject of research <cit.>.
  
  In this paper, we study the upwind margin of the White Sands Dune Field, where the dune instability leads
  to spatially amplifying sand waves developing downstream <cit.>. We follow the spatiotemporal evolution of incipient dunes and extract their wavelength and propagation velocity, as well as their temporal and spatial growth rates. We then show that these four quantities all quantitatively compare to the predictions of spatial linear stability analysis.
  
  
  
  
  
  
  
  
  
  § WHITE SANDS DUNE FIELD
  
  
  
  White Sands Dune Field is located in southern New Mexico, USA. The sand covers an area of about 400 km^2, resulting in the largest gypsum dune field on Earth. Dominant winds are mainly toward the northeast and shape the sedimentary bed into transverse and barchan dunes, progressively turning into parabolic dunes as the vegetation cover increases (Figures <ref>a and <ref>b) <cit.>. Dunes emerge on the upwind margin (Figure <ref>c). There, the sediment is made of coarse, elongate, and angular grains (see Figure S2 in the supporting information), with a measured diameter d = 670 ± 120 μm and a bulk density ρ_ p = 2300 ± 100 kg m^-3. The saturation length, relevant in the process of dune emergence (see section <ref>), directly depends on these grain properties <cit.>:
  
  
  
      L_ sat≃  2.2  ρ_ p/ρ_ f d = 2.8 ± 0.5  m,
  
  
  
  where ρ_ f≃ 1.2 kg m^-3 is the air mass density in ambient conditions. As one moves further into the dune field, the grain diameter and angularity both decrease, due to abrasion and aeolian sorting <cit.>. Because we restrict our analysis of the dune development to the first kilometer along the margin, the grain characteristics are assumed to be homogeneous. The grain roughness leads to a measured avalanche slope μ = 0.8 ± 0.05 (see supporting information section S1).
  
  The sand flux is calculated from the hourly wind data of the weather station at Holloman Air Base (KHMN, 32^∘51'N, 106^∘06'W), using the method described in Cour14 (see also supporting information section S2). The wind is characterized by its shear velocity u_*, representative of the logarithmic profile inside the turbulent boundary layer. Its threshold value u_ th below which saltation cannot sustain steady transport is estimated with u_ th = a √((ρ_ p/ρ_ f) gd) = 0.35 m s^-1, and a ≃ 0.1 <cit.>. The corresponding typical velocity ratio u_*/u_ th is then about 1.26 ± 0.05 (all these values are gathered in Table S1). As shown in Figure <ref>a and documented by Pede15, the wind regime is multimodal. Southwesterly winds dominate the transport as noted by the nearly unimodal sand flux distribution toward the northeast. The other modes from the north and southeast moderately impact the dune shape and migration<cit.>.
  
  Elevation data of the blue area in Figure <ref>a have been obtained using lidar-derived topography at five different times (June 2007, June 2008, January 2009, June 2010 and August 2015). Along the upwind margin, we extracted 75 dune profiles from the surface elevation data, spaced 50-m apart and aligned with the direction of the resulting sand flux (Figure <ref>b). The average topography is removed using a Butterworth high-pass filter (order = 5, cutoff frequency = 0.005 m^-1). We limit our analysis to the first dunes of the filtered profile to ensure that we stay as much as possible in the early stage of dune growth (red area in Figure <ref>c). As shown by figure <ref>d, these incipient dunes have very low aspect ratios and do not exhibit any slipfaces <cit.>. The detrended bed elevation exhibits a spatially amplified oscillating behaviour (Figure <ref>a and <ref>b). We now interpret these profiles using the theoretical framework provided by the linear stability analysis developed in the next section.
  
  
  
  
  
  
  
  
  
  
  § DUNE LINEAR STABILITY ANALYSIS
  
  
  
  Here, we consider a unidirectional wind blowing at a constant shear velocity u_*, over a flat sedimentary bed. Above the transport threshold velocity u_ th, the saturated sediment flux induced by this flow follows a quadratic law:
  
  
  
      q_ sat = Ω(u_ *^2 - u_ th^2),
  
  
  
  where Ω is a dimensional constant that depends on fluid and grain properties <cit.>. In natural conditions, however, a sedimentary bed is never perfectly flat nor infinite, and these irregularities can be seen as the sum of different perturbations. The purpose of the linear stability analysis is precisely to study the temporal or spatial evolution of the bed in response to a perturbation of a given time or length scale. The emerging dune pattern is then expected to be dominated by the most unstable scale, associated with a sinusoidal mode of wavelength λ and propagation velocity c, and whose amplitude grows in time with a rate σ or in space over a length Λ.
  
  Above such a sinusoidal bed, wind and sand flux are also modulated. As described in the introduction, the basal shear stress is not in phase with the topography; this is quantified with two dimensionless coefficients, 𝒜 and ℬ, which represent the in-phase and in-quadrature components, respectively. They are weak (logarithmic) functions of λ <cit.> but can be, in the first approximation for our purpose, considered as constant parameters of the model. The corresponding upwind shift of the wind with respect to the bed elevation is ∼λℬ/(2π𝒜). Similarly, the actual sediment flux q is not saturated but delayed with respect to the basal shear stress, a process quantified by the saturation length L_ sat. These are the main physical mechanisms of the dune formation model from which the stability analysis is derived (see supporting information section S3 for the proper technical derivation and related theoretical figures).
  
  
  Consider first a large flat extent of sand. Under the action of the wind, dunes emerge everywhere simultaneously: there is no spatial development of the pattern. A spatial sinusoidal perturbation characterized by a given wavelength λ and an initial amplitude C_0 can grow or decay in time, at a rate σ, in response to the wind shift and the flux lag. The perturbation also propagates at a velocity c. The bed elevation along the direction x of the wind can be written as:
  
  
  
      h(x,t) = C_0  e^σ tcos[2π/λ(x-c t)].
  
  
  
  Both temporal growth rate and propagation speed can be computed as functions of λ from the analysis of the equations coupling the flow, sediment transport and bed evolution, constituting the dispersion relation of sand waves (supporting information section S3). Positive values of the growth rate are associated with unstable perturbations, and these are typically for large values of λ. Conversely, small wavelengths are stable, with σ<0.
  
  
  Consider now a sediment bed bounded upwind such that disturbances cannot grow at a specific position in space, noted here x=0. Dunes emerge from the selective amplification of disturbances propagating downwind from the field entrance. This results in a spatial development of the pattern. There is no temporal growth: at a given location, the amplitude of the bed oscillation is the same at any time. The form of a sinusoidal mode of initial amplitude C_0 can be written as:
  
  
  
      h(x,t) = C_0  e^x/Λcos[2π/λ(x-c t)],
  
  
  
  where Λ^-1 is the spatial growth rate of the dunes.
  
  Neutral modes are the same in both spatial and temporal analysis. Their wavelength λ_c is characterised by vanishing growth rates σ(λ_c)=0 and Λ^-1(λ_c)=0, such that,
  
  
  
      λ_c = 2π𝒜/ℬ - 1/μ(u_ th/u_*)^2  L_ sat,
  
  
  
  and separates growing (λ > λ_c) from decaying (λ < λ_c) perturbations. It can thus be interpreted as a minimal dune size.
  
  Performing the temporal linear stability analysis (denoted by subscript T), in the limit L_ sat/λ_c ≪ 1, the characteristics of the fastest growing perturbation read:
  
  
  
      λ_ T   ∼   3/2λ_c,
          
      σ_ T   ∼   Q/L_ sat^2 𝒜/2 (2 π L_ sat/λ_ T)^3,
          
      
      	c_ T   ∼   Q/L_ sat 𝒜 2 π L_ sat/λ_ T,
  
  
  
  where Q = Ω u_ *^2 gives the characteristic scale of the sediment flux associated with the wind regime <cit.>.
  
  Conversely, spatial growth rate reaches a maximum 1/Λ_S at some specific value of the wavelength, noted λ_ S, corresponding to the propagation velocity c_ S. Unfortunately, no simple analytical and accurate formulae like (<ref>-<ref>) can be derived for these quantities (see supporting information section S3). Temporal and spatial analyses are consistent, and we typically find λ_ S≃ 1.3 λ_ T and c_ S≃ 0.77 c_ T. The numerical factors in these relations do not vary by more than a few percent upon changing the parameters 𝒜 and ℬ within a reasonable range of values. In the spatial development of the instability, an individual bump grows in height while propagating downwind at a constant velocity. Its amplitude therefore varies exponentially with respect to time, and one can define a temporal growth rate ς_ S≡ c_ S/Λ_ S (inset of Figure <ref>c) in the frame of reference of the bump i.e. a Lagrangian growth rate. The theoretical analysis provides the approximate relation ς_ S≃ 0.43σ_ T. Both ς_ S and Λ_S can be measured from the field data, providing equivalent information.
  
  
  
  
  § FIELD DATA ANALYSIS
  
  
  After removing the average topography (Figure <ref>c), the detrended bed elevation profile exhibits an exponentially amplifying sinusoidal shape, as predicted by the spatial linear theory for dune emergence (Figure <ref>a). Using these profiles and their temporal evolution (Figure <ref>b), we have extracted the three independent characteristics of the pattern (λ, c and Λ or ς) using two different methods. For each profile, we either look at each peak separately (peak-to-peak method), or extract quantities averaged over the whole profile (global approach).
  
  The wavelength λ is computed by autocorrelation of the bed elevation profile (global method), and from the spacing between two adjacent peaks (peak-to-peak method). The fit of an exponential to the peaks of each profile gives the spatial growth length Λ (global method, see dashed line in Figure <ref>a). The spatial growth length is alternatively computed from the difference in height between two adjacent peaks (peak-to-peak method, see Figure <ref>c).
  
  The Lagrangian growth rate ς and propagation velocity c are obtained by fitting exponential and linear functions to the temporal variation of the dune height and position, respectively (see equations 43 and 44 of the supporting information section S3). The peak-to-peak method looks at the height and position of each peak separately (see inset of Figure <ref>c). For a global measurement, the average propagation speed can be determined from the cross-correlation curve between the same profiles at different times, and the average Lagrangian growth rate from the temporal evolution of the bed elevation standard deviation.
  
  Importantly, to extract the average values of ς and c, we take into account the temporal variations of the characteristic sand flux Q. Indeed, both quantities vary in time proportionally to Q.  We also remove periods of time when the wind is below the transport threshold (see supporting information section S4).
  
  
  
  
  
  
  § TIME AND LENGTH SCALES OF THE INCIPIENT DUNES
  
  
  
  The output of the analysis of the 75 transects is shown in Figure <ref>. Both peak-to-peak and global methods exhibit similar distributions for the wavelength λ, the propagation velocity c, the Lagrangian growth rate ς and the growth length Λ, with clear dominant (most probable) values. The incipient dune wavelength and growth length are both on the order of a hundred meters; their propagation velocity is around 5 m yr^-1 and their growth rate is about 0.015 yr^-1. These values, as well as the typical dispersion around them (i.e. the width of these distributions), are more precisely reported in Table S2. For λ and c, our results are consistent with the measurements of Phil19, made on a single elevation profile. Their dispersion is on the order of 20%, because these quantities can be measured with a good accuracy, especially with the global method using correlation. As the measurement of Λ and ς is more delicate, the corresponding distributions are more dispersed. The peak-to-peak method is actually sensitive to the behaviour of individual peaks, that can respond to various types of local disturbances. For example, they may induce irregularities in the spacing between the peaks, or asymmetry between positive and negative detrended topography (Figure <ref>b). As a result, a few negative values of the velocity, growth rate and characteristic growth length are reported. Nevertheless, these data provide reliable and meaningful statistics to test the theory, which must be able to account for those four quantities concomitantly. The free parameters of the linear analysis are the hydrodynamic coefficients 𝒜 and ℬ, as the others are set independently with the wind and sediment properties (see Table S1).
  
  
  The incipient dune wavelength peaks around 120 m, and is therefore significantly larger than the usually reported value (≃ 20 m, for sand particles of size 180 μm) <cit.>. How can the theory reproduce such a large value? First, as the wavelength is proportional to the saturation length, and thus to the grain size (equations <ref>, <ref> and <ref>), the presence here of much coarser grains provides a factor 670/180 ≃ 3.7 corresponding to the ratio of grain diameter. Second, the most unstable wavelength is predicted to increase when sediment transport occurs close to the threshold, in relation to the denominator of equation <ref> <cit.>. Here, with u_*/u_ th≃ 1.3, we can expect an additional factor of 2 with respect to a situation far above threshold. The simultaneous fit of the four quantities predicted by the spatial linear stability analysis to the data allows us to reproduce quantitatively the dominant wavelength and growth rate. The predicted values fall in the peak of the distributions (red lines in Figures <ref>a and <ref>c). This adjustment however overestimates the growth length and the propagation velocity, whose predictions exceed the dominant values by an amount comparable to the peaks' width (red lines in Figures <ref>b and <ref>d).
  
  This discrepancy can be understood by questioning our approximation of a unidirectional wind. A finer analysis of the flux rose shows in fact secondary winds, with non-zero components perpendicular to the crest toward the southwest. Reversing winds have cumulative effects on the growth rate and the selected wavelength. They however partially cancel each other out, impeding the propagation, and thus the spatial development of the dune pattern. Such a process is supported by observations of reactivation surfaces formed by reversing winds in the stratigraphy at White Sands <cit.>. The value of the characteristic sand flux Q (given in Table S1) is a time average that does not account for changes in wind orientation. The ratio between the scalar and vector averages of Q, taking into account the variations in orientation of the sand fluxes over time, is about 0.6 (see supporting information section S4). Once corrected by this ratio, the predicted migration velocity and growth length come into quantitative agreement with the corresponding dominant values of the distributions (orange lines in Figure <ref>).
  
  
  
  
  
  
  
  
  
  
  § DISCUSSION
  
  
  
  The exponentially amplified sinusoidal behaviour of the White Sands bed elevation profiles is field evidence for the spatial development of the dune instability. The linear analysis is able to quantitatively reproduce the three characteristics of the emergent pattern: dune wavelength, propagation velocity and growth length (or equivalently, growth rate). To obtain this agreement, the two hydrodynamic coefficients where adjusted, resulting in 𝒜 = 3.6 ± 0.6 and ℬ = 1.9 ± 0.3. All the other parameters of the theory were fixed independently from sediment and wind data, either by direct measurement (grain diameter and density, avalanche slope) or using well calibrated relationships (saturation length, sediment flux). Note that the uncertainty on the determined coefficients is dominated by the those of the fixed parameters, rather than by the dispersion in the measurements of c, λ and Λ. Overall, the value of L_ sat, involved in λ and Λ, directly affects the estimate of ℬ/𝒜, while the value of Q, involved in c, mostly affects 𝒜.
  
  Importantly, the concomitant agreement of c, λ and Λ is a stringent test of the theory. This study is therefore a step forward in the general `dune inverse problem', trying to infer, for example, grain or wind properties from dune characteristics <cit.>. It is remarkable that the resulting values determined from our field data are very close to those directly measured by Clau13 on a single 40 m long dune (𝒜≃ 3.4 and ℬ≃ 1.6), as well as to the predictions of hydrodynamic models <cit.>. This study thus confirms the reliability of the linear analysis in the interpretation of field, experimental or numerical data relating to the emergence of incipient dunes <cit.>.
  
  A limitation of the linear theory is of course the presence of non-linear effects. They occur when the aspect ratio of the sand waves becomes too large or when the dunes interact with each other, so that each bed perturbation cannot be considered as independent of the others. Bumps with aspect ratios above ≃ 1/13 are expected to start to develop flow recirculation on their downwind side, usually associated with the formation of an avalanche slip face <cit.>. In the region we have studied, the waves furthest inside the dune field could reach aspect ratios of about 1/10, but no slip faces were observed. Similar to Phil19, we also recognize the coexistence of multiple wavelengths at the upwind side of the profiles (associated with different celerities and growth rates or lengths), and these are partly the cause of the distribution widths in Figure <ref>. We could not, however, infer from these data signs of interactions, such as collisions, coalescence or ejection <cit.>.
  
  Although studied here at the boundary of a dune field, the spatial development of the dune instability is also present on pre-existing large dunes, providing a similar upwind boundary condition in terms of sand availability. As a matter of fact, the emergence of bed oscillations on the flanks of barchans has been proposed as a key mechanism to understand their stability, as these superimposed waves eventually grow until they can break from the horns, causing large sand losses <cit.>. Likewise, in narrow bidirectional wind regimes, the growth of the instability over elongating linear dunes breaks them into trains of barchans <cit.>. This work therefore provides a reliable base to study the stability of large dunes and thus the formation of large-scale structures inside dune fields <cit.>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  This collaborative work was initiated at the International Conference on Aeolian Research
  (ICAR), 2018, Bordeaux, France. We acknowledge financial support from the Univ-EarthS LabEx program of Sorbonne Paris Cité (grants ANR-10-LABX-0023 and ANR-11-IDEX-0005-02) and the French National Research Agency (grant ANR-17-CE01-0014/SONO). Clément Narteau acknowledges support from the National Science Center of Poland (grant 2016/23/B/ST10/01700). Philippe Claudin acknowledges support from TOAD (The Origine of Aeolian Dunes) project as an external partner (Natural Environment Research Council, UK and National Science Foundation, USA; NE/R010196NSFGEO-NERC). Ryan C. Ewing acknowledges support from the White Sands National Monument through NPS-GC-CESU Cooperative Agreement #P12AC51051.
  Meteorological data used in this manuscript are hosted by MesoWest under the station code `KHNM' (<https://mesowest.utah.edu/cgi-bin/droman/mesomap.cgi?state=NM   rawsflag=3>).
  Topographic data can be found in the public repositories Texas Data Repository and OpenTopography (2007-doi:10.18738/T8/WUNF0G, 2008-doi:10.18738/T8/HQVUSX, 2009-doi:10.5069/G9Q23X5P, 2010-doi:10.5069/G97D2S2D, 2015:<https://portal.opentopography.org/usgsDataset?dsid=USGS_LPC_NM_WhiteSands_2015_LAS_2017>).
  Supplementary figures, tables, and text can be found in the supporting information. We thank the anonymous reviewers for their careful reading of our manuscript and their insightful comments and suggestions.
  
  
  
  
  
  
  
  
  
  
  
  
  