# Modelling cargo transport in crowded environments: effect of motor
association to cargos

Sutapa Mukherji<EMAIL_ADDRESS>Dhruvi K. Patel  Mathematical
and Physical Sciences division, School of Arts and Sciences, Ahmedabad
University, Navrangpura, Ahmedabad 380009, India

###### Abstract

In intracellular transports, motor proteins transport macromolecules as cargos
to desired locations by moving on biopolymers such as microtubules. Recent
experiments suggest that cargos that can associate motor proteins during their
translocation have larger run-length, association time and can overcome the
motor traffic on microtubule tracks. Here, we model the dynamics of a cargo
that can associate at the most $m$ free motors present on the track as
obstacles to its motion. The proposed models display competing effects of
association and crowding, leading to a peak in the run-length with the free
motor density. This result is consistent with past experimental observations.
For $m=2$ and $3$, we show that this feature is governed by the largest
eigenvalue of the transition matrix describing the cargo dynamics. In all the
above cases, free motors are assumed to be present as stalled obstacles. We
finally compare simulation results for the run-length for general scenarios
where the free motors undergo processive motion in addition to binding and
unbinding to or from the microtubule.

##  1 Introduction

Intracellular transport often involves directional movements of motor proteins
on biopolymers such as microtubules or actin filaments howard1 ; schliwa .
Three major classes of motor proteins known as kinesin, dynein and myosin are
responsible for such transports. Using the energy derived from the hydrolysis
of adenosinetriphosphate (ATP) molecules, motor proteins transport different
types of cargos such as cellular organelles, protein complexes, mRNAs etc. to
desired locations in the cell. Such cargo movements are essential for various
cellular functions such as cell morphogenesis, cell division, cell growth etc.
This motion is processive in the sense that motor proteins typically move over
several successive steps before detaching from the microtubule. Early studies
block ; cross ; howard2  on intracellular transport revealed the underlying
mechanism behind motor transport and how various properties such as the run-
length, velocity etc. depend, for example, on the external force or the
concentration of ATP molecules. While many of these studies are around the
transport by a single motor, it is believed that cargos are often transported
by multiple motors gross ; holzwarth ; gelfand ; unger  which help cargos
remain bound to the biopolymer for a longer time. Experimental and theoretical
studies holzwarth ; klumpp ; beeg  show that the presence of several motors
helps the cargo overcome the viscous drag of the cytoplasm and have larger
velocity as compared to transport by single motors. The cooperation of several
motors also leads to longer run-length of the cargo before it detaches from
the microtubule. Further, in vitro experiments indicate that transport
processes by multiple motors can be efficiently regulated by controlling the
number of engaging motors vershinin . Besides these studies, there have been
extensive experimental and theoretical studies attempting to understand the
collective nature of transports involving many motors under diverse conditions
leduc ; segregation ; clogging ; helical ; tasep1 ; multi-opposite ; frey-lang
; multi-same ; laneswitch ; gov ; freygait .

Quite often such transport processes take place in a crowded environment of
the intracellular space. This is in particular true for the axon region of the
neuron cell where a dense network of biopolymers, pre-exisiting organelles and
the narrow geometry of the axon together give rise to a crowded environment
that can impede cargo movements. However, despite crowding, it is found that
the cargo transport happens in a robust manner without significant jamming or
cargo dissociation. Experiments elucidating cargo transport in crowded
environments indicate that motor proteins can adapt alternative strategies
that might help them circumvent the crowding problem leduc ; surrey . In a
recent experimental study aimed at understanding the motion of a cargo in a
crowded environment, Conway et al conway ; conway1  studied the motile
properties of quantum dot (Qdot) cargos, that can associate multiple kinesins,
on a microtubule crowded with free kinesin motors. While comparing the motile
properties of free kinesins and the Qdot cargos in crowded conditions, cargos
were found to display longer run-lengths and association times as compared to
free kinesins as the motor density increased. This difference prompted the
prediction that the property of a cargo to associate multiple motors helps
increase its run-length, association time and overcome the motor traffic. It
was observed that while translocating, Qdot cargos could associate kinesins
from the microtubule pool, dissociate kinesins attached to itself, or
associate kinesins that are already moving along the microtubule and move
together subsequently.

Motivated by this work, here, we propose mathematical and computational models
to characterise the motion of a single cargo on a track crowded by free
kinesin motors. During its translocation, the cargo can associate free motors
which impede the motion of the cargo by occupying the forward sites on the
microtubule. We assume that upon such association, a kinesin detaches itself
from the microtubule rendering the forward site free. Our aim is to find how
the interplay of the kinesin-association property of the cargo and the
crowding along the track affects cargo’s motile properties, for example, its
run-length and association time etc. To our knowledge, this is the first
modelling study of cargo transport where the cargo has the ability to
associate kinesins present on a crowded microtubule track.

To this end, we study cargo transport under different scenarios described
below. (1) In the simplest scenario, we assume that the cargo is always bound
to the microtubule. Along the path of the cargo, the microtubule binding sites
are randomly occupied by free kinesin motors. We assume that the kinesin that
gets associated to the cargo during its translocation plays no specific role
in facilitating the forward motion of the cargo other than freeing the forward
site. This is equivalent to assuming that the cargo removes the kinesin
molecule occupying the forward site via the association process. For this
model, referred as "Model 1" below, we find the average velocity of the cargo.
(2) In model 2, we assume that the cargo can bind more than one, say, at the
most $m$ number of kinesins. This is based on the predictions that the cargo
may have a finite number of kinesin binding sites conway . Hence, the cargo
can associate a kinesin occupying the forward site provided it has a free
binding site available. We consider $m=2,\ 3,{\rm and}\ 4$ in the following
analysis. Kinesins attached to the cargo can detach from it and a free kinesin
from the intracellular space can attach to the cargo at given rates. Finally,
we implement the condition that a cargo can no longer be on the microtubule
track if all the kinesin molecules detach from the cargo. A generalised
version of the mathematical formulation of model 1 allows us to analyse the
cargo motion obeying above rules for $m=2$ and $3$. Finally, run-lengths of
$m=2,\ 3,{\rm and\ 4}$ are found upon numerically simulating the cargo
dynamics. The motion of the cargo following different dynamical rules are
shown in figure 1.

Figure 1: (A) Cargo transport on a microtubule in the presence of free
kinesins bound to microtubules. The microtubule is represented by a one-
dimensional lattice with a site representing a tubulin dimer, the basic
subunit of a microtubule. (B) A cartoon of a cargo bound to two kinesins. The
cargo has more than one kinesin binding sites. (C) The process of kinesin
association to cargo at rate $r_{an}$. (D) The process of attachment of a
kinesin from the intracellular space to the cargo at rate $\omega_{a}$. (E)
The process of kinesin detachment from the cargo at rate $\omega_{d}$.

(3) In models 1 and 2, free kinesins are assumed to be stalled on the
microtubule. In model 3, we simulate cargo dynamics in the presence of moving
kinesins as well as random processes of kinesin binding and unbinding to or
from the microtubule. We compare the run-length of the cargo (with $m=3$) in
the presence or absence of various processes mentioned above.

##  2 Models and Results

###  2.1 Model 1

The motion of the cargo is modelled considering the following dynamical rules.
(a) The cargo transported by a kinesin starts its forward journey from a given
point on a one-dimensional track (often referred below as a lattice)
representing the microtubule. (b) For all the lattice sites ahead, we assume
an initial, random distribution of free kinesins. The average kinesin density
on the lattice is represented by $r_{m}$. These kinesins are assumed to be
stalled. (c) The cargo moves to the forward site provided the forward site is
not occupied by a free kinesin. (d) If the forward site in front of the cargo
is occupied by a free kinesin, the cargo can associate the kinesin with itself
at rate $r_{an}$ rendering the forward site free.

In order to build the mathematical model, we consider possible configurations
that two neighbouring sites can have when the first one of them is occupied by
the cargo. For $i$ and $(i+1)$-th sites, with the cargo being at the $i$-th
site, the $(i+1)$-th site can be either empty or occupied by a free kinesin
molecule. We denote the probability of finding $(i+1)$-th site empty with the
$i$-th site occupied by the cargo at time $t$ by $P(i,t)$. Similarly, the
probability of finding $(i+1)$-th site occupied by a free kinesin while the
cargo is at the $i$-th site at time $t$ is $Q(i,t)$. Figure (2) shows these
configurations as well as possible transitions from one configuration to the
other as the cargo translocates forward. The following equations describe how
these two configurations evolve with time sm .

$\displaystyle\frac{dP(i)}{dt}=r_{an}Q(i)+(1-r_{m})P(i-1)-P(i),$ (1)
$\displaystyle\frac{dQ(i)}{dt}=-r_{an}Q(i)+r_{m}P(i-1).$ (2)

The first term on the RHS of equation (1) indicates a cargo-association
process due to which a Q-type configuration transitions to a $P$-type
configuration. The term with the pre-factor $(1-r_{m})$ indicates the motion
of the cargo from $(i-1)$-th site to $i$-th site while the $(i+1)$-th site is
vacant. While the forward motion happens with unit rate, the factor
$(1-r_{m})$ indicates the probability that after the forward motion,
$(i-1)\rightarrow i$, the cargo lands in a P-type configuration i.e. the
$(i+1)$-th site is unoccupied by a kinesin. The last term in (1) is a loss
term which indicates that a cargo has moved from the $i$-th site to the
$(i+1)$-th site. In equation (2), the first term on the RHS indicates a loss
of a $Q$-type configuration due to the association process. The second term is
a gain term due to the hopping of the cargo from $(i-1)\rightarrow i$ where
$(i+1)$-th site is occupied by a free kinesin.

Figure 2: (A) Two possible configurations of two neighbouring sites with the
first site occupied by the motor transporting the cargo. (B) Possible
transitions associated with the transport of the cargo to the next site.

To find the average properties of the cargo motion, we define generating
functions corresponding to the two probabilities as

$\displaystyle\tilde{P}(\gamma)=\sum_{i=-\infty}^{\infty}\gamma^{i}P(i),\ {\rm
and}\ \tilde{Q}(\gamma)=\sum_{i=-\infty}^{\infty}\gamma^{i}Q(i).$ (3)

In terms of these generating functions, the time evolution equations are

$\displaystyle\frac{d}{dt}\tilde{P}(\gamma)=r_{an}\tilde{Q}(\gamma)+(1-r_{m})\gamma\tilde{P}(\gamma)-\tilde{P}(\gamma),\
{\rm and}$ (4)
$\displaystyle\frac{d}{dt}\tilde{Q}(\gamma)=-r_{an}\tilde{Q}(\gamma)+r_{m}\gamma\tilde{P}(\gamma).$
(5)

The average position of the cargo can be found from the probabilities as

$\displaystyle\langle
i\rangle=\sum_{i=-\infty}^{\infty}i[P(i)+Q(i)]=\frac{d}{d\gamma}[\tilde{P}(\gamma)+\tilde{Q}(\gamma)]\mid_{\gamma=1}.$
(6)

The average velocity of the cargo is obtained from $v=\langle i\rangle/t$
where $t$ is the time taken to travel an average distance $\langle i\rangle$.
Solving equations (4) and (5), the average velocity of the cargo is found as
(see Appendix A for details)

$\displaystyle v=\frac{r_{an}}{r_{an}+r_{m}}.$ (7)

###  2.2 Model 2

Here we generalize the mathematical framework, discussed in the previous
section, for higher values of $m$ taking into account the possibilities of
detachments of the cargo from the microtubule. In the following, we discuss
the mathematical model for $m=2$ and simulation results for $m=2,\ 3,\ {\rm
and}\ 4$. The cargo dynamics for $m=3$ is discussed in Appendix B.

For $m=2$, the cargo has two kinesin binding sites. Hence it can associate at
the most two kinesins. The basic rules for cargo transport in this case are
listed below. (a) As before, we begin with an initial, random distribution of
stalled free kinesins on a one-dimensional lattice. The average density of
free kinesins is $r_{m}$. (b) The cargo attached to a kinesin starts its
forward journey from a given point on the lattice. (c) If the forward site is
blocked by a free kinesin, the cargo can associate the kinesin with itself at
a rate $r_{an}$ provided the cargo has only one kinesin bound to it. (d) A
kinesin bound to the cargo can detach from the cargo at rate $\omega_{d}$ and
a free kinesin from the intracellular space can bind to the cargo at rate
$\omega_{a}$ provided the cargo has only one kinesin attached to it. (e) A
cargo is not attached to the microtubule if all the kinesins detach from the
cargo. Thus, we are not being specific about how many kinesins are actively
transporting the cargo or how many remain bound to the cargo without
participating in cargo transport actively.

As before, we begin with two possible configurations of, say, $\\{i,i+1\\}$-th
sites where $i$-th site is occupied by the cargo. However, here the cargo can
be in two possible states - bound to one kinesin or bound to two kinesins.
Hence, the probabilities are defined in the following way. $P_{n}(i,t)$
($n=1,\ 2)$ represents the probability, at time $t$, of the configuration
where the cargo, located at $i$-th site, is bound to $n$ kinesins and the
$(i+1)$-th site is empty. Similarly, $Q_{n}(i,t)$ ($n=1,\ 2$) represents the
probability, at time $t$, of the configuration where the cargo, located at
$i$-th site, is bound to $n$ kinesins and the $(i+1)$-th site is occupied by a
free kinesin.

The probabilities of various configurations change with time as per the
following equations,

$\displaystyle\frac{d}{dt}P_{2}(i)=(1-r_{m})P_{2}(i-1)-P_{2}(i)+r_{an}Q_{1}(i)-\omega_{d}P_{2}(i)+\omega_{a}P_{1}(i),$
(8)
$\displaystyle\frac{d}{dt}P_{1}(i)=(1-r_{m})P_{1}(i-1)-P_{1}(i)-\omega_{a}P_{1}(i)+\omega_{d}(P_{2}(i)-P_{1}(i)),$
(9)
$\displaystyle\frac{d}{dt}Q_{2}(i)=r_{m}P_{2}(i-1)+\omega_{a}Q_{1}(i)-\omega_{d}Q_{2}(i),\
{\rm and}$ (10)
$\displaystyle\frac{d}{dt}Q_{1}(i)=r_{m}P_{1}(i-1)-r_{an}Q_{1}(i)+\omega_{d}(Q_{2}(i)-Q_{1}(i))-\omega_{a}Q_{1}(i).$
(11)

The $r_{\rm an}$-dependent term in equation (8) represents a process of
kinesin association by the cargo. Due to this process, a $Q_{1}$-type
configuration transitions to a $P_{2}$-type configuration.
$\omega_{a}$($\omega_{d}$) dependent terms represent kinesin
attachment(detachment) processes to(from) the cargo. For example, the
$\omega_{d}$ dependent term in equation (8) represents detachment of a kinesin
due to which the cargo transitions from $P_{2}$ state to $P_{1}$ state. In
addition to above equations, we introduce probabilities $P_{0}(i,t)$ and
$Q_{0}(i,t)$ of having situations where the cargo bound to one kinesin
residing at $i$-th site at time $t$ loses its kinesin. These probabilities
change with time as per the equations

$\displaystyle\frac{d}{dt}P_{0}(i)=\omega_{d}P_{1}(i)\ {\rm and}\
\frac{d}{dt}Q_{0}(i)=\omega_{d}Q_{1}(i).$ (12)

Defining generating functions as
$\tilde{P}_{n}(\gamma,t)=\sum_{i=-\infty}^{\infty}\gamma^{i}P_{n}(i,t)$, and
$\tilde{Q}_{n}(\gamma,t)=\sum_{i=-\infty}^{\infty}\gamma^{i}Q_{n}(i,t)$ (where
$n=0,\ 1,\ 2$), we can rewrite equations (8)-(11) as

$\displaystyle\frac{d}{dt}{\bf H}(\gamma,t)={\bf S}{\bf H}(\gamma,t),$ (13)

where ${\bf H}$ is a column matrix

$\displaystyle{\bf H}(\gamma,t)={\begin{pmatrix}&\tilde{P}_{2}(\gamma,t)\\\
&\tilde{P}_{1}(\gamma,t)\\\ &\tilde{Q}_{2}(\gamma,t)\\\
&\tilde{Q}_{1}(\gamma,t)\\\ \end{pmatrix}}$ (14)

and ${\bf S}$ is a $4\times 4$ matrix

$\displaystyle{\bf
S}={\begin{pmatrix}&(1-r_{m})\gamma-1-\omega_{d}&\omega_{a}\ &0&r_{an}\\\
&\omega_{d}&(1-r_{m})\gamma-1-\omega_{a}-\omega_{d}&0&0\\\
&r_{m}\gamma&0&-\omega_{d}&\omega_{a}\\\
&0&r_{m}\gamma&\omega_{d}&-(\omega_{a}+\omega_{d}+r_{an})\\\ \end{pmatrix}}.$
(15)

In order to have an estimate of the association time of the cargo and how it
is impacted by various processes, we have studied the quantity
$[\tilde{P}_{0}(\gamma,t)+\tilde{Q}_{0}(\gamma,t)]\mid_{\gamma=1}$. This
quantity being identical to $\sum_{i=-\infty}^{\infty}[P_{0}(i,t)+Q_{0}(i,t)]$
indicates the total probability of cargo being left with no kinesin bound to
it while being at any point on the lattice. A plot of this quantity for
different parameter values are shown in figures (3) and (4). Over large time,
this quantity approaches unity indicating cargo losing all the kinesins
leading to the detachment of the cargo from the microtubule. The approach of
this quantity to unity is what provides us with an estimate of the association
time of the cargo to the microtubule. A fast approach to unity indicates a low
association time of the cargo. For both figures, we have chosen the same sets
of values for the kinesin-association rate, $r_{an}$. The increase or decrease
in $\omega_{a}$ and $\omega_{d}$, respectively, is expected to increase the
association time of the cargo. Figures show that reducing the kinesin
detachment rate, $\omega_{d}$, from the cargo has much stronger effects on the
association time as compared to increasing the kinesin attachment rate,
$\omega_{a}$.

Figure 3: $y-$axis represents the total probability
($\sum_{i=-\infty}^{\infty}[P_{0}(i,t)+Q_{0}(i,t)]$) of the cargo losing its
last kinesin at time $t$ while being at any site on the lattice. For this
plot, $r_{m}=0.5$ and $\omega_{a}=0.01$. Figure 4: $y-$axis represents the
total probability ($\sum_{i=-\infty}^{\infty}[P_{0}(i,t)+Q_{0}(i,t)]$) of the
cargo losing its last kinesin at time $t$ while being at any site on the
lattice. For this plot, $r_{m}=0.5$ and $\omega_{d}=0.01$.

In figure (5), we have shown how the total probability of cargo detachment at
any point on the lattice is influenced by the kinesin density, $r_{m}$. The
figure shows that at a low value of the kinesin-association rate by the cargo,
$r_{\rm an}$, the extent of crowding influences the cargo-association time to
the microtubule only mildly. The situation changes significantly when the
kinesin-association rate is high. In this case, the association time of the
cargo to the microtubule, in general, increases significantly. Further, for
large $r_{\rm an}$, the crowding density of free kinesins affects the
association time of the cargo significantly with the association time being
larger for lower crowding density.

Figure 5:  $y-$axis represents the total probability
($\sum_{i=-\infty}^{\infty}[P_{0}(i,t)+Q_{0}(i,t)]$) of the cargo losing its
last kinesin at time $t$ while being at any site on the lattice. For this
plot, $\omega_{a}=\omega_{d}=0.001$.

The dependence of the run-length of the cargo on the crowding density can be
obtained upon solving equations (8)-(11) numerically. Figure (6) shows run-
length plots for different values of the kinesin-association rate, $r_{\rm
an}$, and kinesin attachment and detachment rates, $\omega_{a}$ and
$\omega_{d}$, respectively. For small $r_{\rm an}$, the run-length decreases
monotonically. However, for large $r_{\rm an}$, the run-length increases
initially for low crowding. In this case, due to large $r_{\rm an}$, the cargo
benefits from the kinesin-association process at low crowding. As the crowding
density increases, due to limited number of binding sites, the cargo no longer
benefits from kinesin association and the run-length decreases. This variation
of the run-length with the crowding density is consistent with earlier
experimental predictions conway . With the increase in the kinesin detachment
rate, $\omega_{d}$, the run-length of the cargo decreases significantly.
However, as found earlier, a decrease in the rate of kinesin attachment,
$\omega_{a}$, to the cargo has mild effect on cargo’s run-length.

Figure 6: Run-length of the cargo plotted with crowding density, $r_{m}$. The
cargo can bind at the most two kinesins ($m=2$).

The fact that the run-length of the cargo for high kinesin-association rates
increases with the crowding density initially is consistent with the estimates
obtained from the analysis of the largest eigenvalue of the transition matrix
$S$ and the association time. In the limit of large time, the solutions for
the probabilities are given by

$\displaystyle{\bf H}\approx c_{1}e^{\lambda_{l}t}{\bf X},$ (16)

where $c_{1}$ is a constant, $\lambda_{l}$ is the largest of the four
eigenvalues of the transition matrix ${\bf S}$ with all of them being negative
and ${\bf X}=(x_{1},\ x_{2},\ x_{3},\ x_{4})^{T}$ is the corresponding
eigenvector. The average distance travelled by the cargo and its average
velocity can be obtained from $\langle
i\rangle=\sum_{i=-\infty}^{\infty}i[P_{1}(i,t)+P_{2}(i,t)+Q_{1}(i,t)+Q_{2}(i,t)]\\\
=\gamma\frac{d}{d\gamma}[\tilde{P}_{1}(\gamma,t)+\tilde{P}_{2}(\gamma,t)+\tilde{Q}_{1}(\gamma,t)+\tilde{Q}_{2}(\gamma,t)]\mid_{\gamma=1}$
and $v={\langle i\rangle}/{t}$, respectively. In the large time limit, the
dominant contribution to the velocity is of the form $v\approx[\gamma
c_{1}e^{\lambda_{l}t}\frac{d\lambda_{l}}{d\gamma}\sum_{i=1}^{4}x_{i}]\mid_{\gamma=1}$.
Using $1/\lambda_{l}$ as an estimate of the association time, $t_{\rm assoc}$,
and finding $\frac{d\lambda_{l}}{d\gamma}\mid_{\gamma=1}$ numerically for
given parameter values, we have plotted $t_{\rm
assoc}\frac{d\lambda_{l}}{d\gamma}\mid_{\gamma=1}$ as a function of the
crowding density, $r_{m}$, in figure (7). Plots in figure (7) display similar
trends as found in figure (6) for the run-length. Although $t_{\rm assoc}v$
gives an estimate of the run-length, the variation in the run-length with the
crowding density as seen in figure (6) essentially arises from $t_{\rm
assoc}\frac{d\lambda_{l}}{d\gamma}\mid_{\gamma=1}$. It can be verified
numerically that the variation in the remaining factors in $v$ is almost
negligible over the entire range of $r_{m}$, $[0:1]$.

Figure 7:  $y$-axis represents the product of two factors that dominate the
nature of the run-length. $t_{\rm assoc}$ is the association time of the cargo
and $\frac{d\lambda_{l}}{d\gamma}$ arises while computing the average velocity
of the cargo (see the text). This plot is for $m=2$ with
$\omega_{a}=\omega_{d}=0.001$.

The dynamical equations for a cargo that can bind at the most three kinesins,
i.e. $m=3$, are shown in Appendix B. The variation of $t_{\rm
assoc}\frac{d\lambda_{l}}{d\gamma}\mid_{\gamma=1}$ with the crowding density
as obtained from the analysis of the largest eigenvalue is shown in figure
(14).

Next we simulate the cargo dynamics with the cargo having $m=2,\ 3,{\rm and}\
4$ kinesin binding sites. Figure (8) shows the change in the run-length of the
cargo with free-kinesin density, $r_{m}$, for $m=2,\ 3,\ {\rm and}\ 4$.
Simulations show an initial increase in the run-length with the free-kinesin-
density for $m=3\ {\rm and}\ 4$; a trend that was shown earlier in figure (6).

Figure 8: Run-length of the cargo as a function of free-kinesin density along
the microtubule track. For this figure, $\omega_{a}=\omega_{d}=0.05$,
$r_{an}=0.4$. The results are obtained upon averaging over $500$ samples for
$m=2$ and $3$, and over $1500$ samples for $m=4$.

###  2.3 Processive motion of free kinesins for $m=3$

In this section, we study the motion of the cargo in the presence of free
kinesins which move processively on the microtubule track. In addition,
kinesins from the intracellular environment can attach to the microtubule and
those walking on the microtubule can leave the microtubule at given rates.

Here we simulate this system using the cellular automaton method. As before,
the microtubule is represented by a one-dimensional lattice. We begin with the
cargo positioned at one end of the lattice. The lattice sites are randomly
occupied by free kinesins with an average density, $r_{m}$. The cargo moves
following the association mechanism mentioned earlier. We assume that the
kinesins move unidirectionally in the same direction as that of the cargo. The
motion of the free kinesins follows the rules of the paradigmatic totally
asymmetric simple exclusion process tasep1 . Accordingly, each kinesin can
walk to the neighbouring site forward provided the target site is not occupied
by another kinesin. The attachment and detachment of kinesins are as per the
Langmuir kinetics considered in frey-lang . A kinesin can attach to a lattice
site at rate $\omega_{a,{\rm kin}}$ provided the site is empty and a free
kinesin can detach from the lattice at rate $\omega_{d,{\rm kin}}$. A kinesin
located at the boundary site can exit from the lattice at rate, $\beta$. We
follow random sequential updating scheme with probability $p$ for cargo update
and $1-p$ for updating the rest of the sites. Depending on the site chosen,
the state of the site (or of the cargo) is updated following the
aforementioned rules. The description of various parameters is provided in
Appendix C.

In figure (9), we have plotted run-lengths for three different scenarios, (i)
free kinesins are stalled (static obstacles) (ii) free kinesins are in motion
and (iii) free kinesins are in motion and they can attach to (or detach from)
the microtubule at rates $\omega_{a,{\rm kin}}$ (or $\omega_{d,{\rm kin}}$).
Plots indicate that in case of processive free kinesins (case (ii)), the run-
length of the cargo peaks at a higher crowding density with a higher maximum
value as compared to the stalled case. The run-length reduces significantly in
case of random attachment and detachment of free kinesins (case (iii)). Figure
(15) in Appendix C shows that the attachment processes lower the run-length
significantly. Although, due to increased effective crowding density, the
cargo remains bound to the microtubule for a large span of time by associating
free kinesins, the crowding restricts the run-length. As a result of this, the
run-length attains its maximum value at a lower value of the crowding density,
$r_{m}$.

Figure 9:  Run-length of the cargo as a function of crowding density, $r_{m}$,
under different conditions. For "processive and attachment/detachment", free
kinesins move processively along the microtubule and can also randomly
attach/detach to/from the microtubule. For "processive" case, free kinesins
only move processively without any attachment/detachment dynamics. In the
"stalled" case, free kinesins do not move and do not attach or detach from/to
the microtubule. For "processive and attachment/detachment" plot,
$\omega_{a,{\rm kin}}=\omega_{d,{\rm kin}}=0.01$. For the rest of the cases,
$\omega_{a,{\rm kin}}=\omega_{d,{\rm kin}}=0$. Other parameter values are
$\omega_{a}=\omega_{d}=0.05$, $r_{\rm an}=0.4$, $p=0.1$, and $\beta=0.6$. The
total number of lattice sites is $2000$. Run-lengths are obtained upon
averaging over $4000$ samples.

Figure (10) shows the variation in the run-length as the association rate
$r_{\rm an}$ is changed. The peaks in the run-lengths are similar to what have
been observed earlier.

Figure 10:  Run-length of the cargo as a function of crowding density,
$r_{m}$, for different values of the kinesin-association rate $r_{\rm an}$.
Free kinesins move processively and they can randomly attach or detach to or
from the microtubule. For this plot $\omega_{a,{\rm kin}}=0.0008,\
\omega_{d,{\rm kin}}=0.0016$ leduc , $p=0.1$, and $\beta=0.6$. The total
number of lattice sites is $2000$. Run-lengths are obtained upon averaging
over $4000$ samples.

Figure (11) shows a comparison of how the association time of the cargo to the
microtubule depends on $r_{m}$ for different values of the kinesin-association
rate, $r_{\rm an}$. The association time is expressed as the total number of
discrete time steps of simulation till the cargo leaves the microtubule. An
increase in $r_{\rm an}$ helps cargo stay attached to the microtubule for a
longer span of time while as per figure (12), the velocity of the cargo
decreases monotonically with the crowding density, $r_{m}$. Further, no
significant variation in the velocity is seen with $r_{\rm an}$. The
association time and the velocity vary with $r_{m}$ in such a manner that
their product exhibits a peak at a specific value of $r_{m}$.

Figure 11:  Association time of the cargo to the microtubule as a function of
crowding density, $r_{m}$, for different values of the kinesin-association
rate ($r_{an}$). Free kinesins move processively and they can randomly attach
or detach to or from the microtubule. For this plot $\omega_{a,{\rm
kin}}=0.0008,\ \omega_{d,{\rm kin}}=0.0016$, $p=0.1$, and $\beta=0.6$. The
total number of lattice sites is $2000$. The association time is expressed as
the total number of discrete time steps of simulation till the cargo leaves
the microtubule. Association times are obtained upon averaging over $4000$
samples. Figure 12:  Velocity of the cargo as a function of crowding density,
$r_{m}$, for different values of the kinesin-association rate ($r_{an}$). Free
kinesins move processively and they can randomly attach or detach to or from
the microtubule. For this plot $\omega_{a,{\rm kin}}=0.0008,\ \omega_{d,{\rm
kin}}=0.0016$, $p=0.1$, and $\beta=0.6$. The total number of lattice sites is
$2000$. For a given $r_{m}$, we have computed the velocity for every sample.
The values presented here are obtained upon averaging of $4000$ samples.

##  3 Summary

The motion of cargos on biopolymeric tracks crowded due to free or cargo-bound
motor proteins is a subject of immense experimental and theoretical
investigations. The central goal of these studies is to understand how cargos
manage to overcome the motor traffic in order to transport necessary materials
in a robust manner. Motivated by some of the experimental observations on
translocation of quantum dot cargos in crowded environments, we have modelled
mathematically and computationally the motion of a cargo that can bind
kinesins present along its trajectory on the microtubule.

In the mathematical modelling, the kinesins on the microtubule track are
assumed to be stalled. Besides taking into account the kinesin-association
property of the cargo, our model incorporates the following dynamical rules.
(i) The cargo has a limited number of kinesin-binding sites as a result of
which it can bind at the most a given number of kinesins, (ii) bound kinesins
can detach from the cargo and kinesins from the intracellular space can bind
to the cargo at certain rates and (iii) the cargo leaves the microtubule if
all the kinesins detach from the cargo. Upon finding the cargo velocity for a
toy model where the cargo never leaves the microtubule and keeps moving
forward by removing obstacles via association, we generalize the mathematical
formulation to take into account the aforementioned aspects of the cargo
dynamics. We show that the two features, namely, the crowding along the
microtubule and the ability of the cargo to associate kinesins have competing
effects on the run-length of the cargo. For low crowding density, as the
crowding density increases, the cargo benefits due to its ability to associate
kinesins. This leads to an increase in the run-length with the crowding
density. However, as the crowding density increases further, due to its
limited number of kinesin binding sites, the cargo does not benefit anymore
through kinesin association. As a consequence, the run-length decreases for
large values of the crowding density. This nature of the run-length has been
predicted earlier from experimental observations. We show that this property
of the run-length is governed by the largest eigenvalue of the transition
matrix describing the dynamics of the cargo. The model can be generalized
further to incorporate other features such as reversals of the cargo,
bidirectional movements of the cargo, pausing of the cargo etc. with the
frequencies of such events depending on the crowding density. The present work
lays a foundation for such studies. Additionally, this analysis may also lead
to testable predictions for cargo’s motile properties once appropriate
parameter values are available.

Next, we have simulated cargo transport with processive motion of free
kinesins as well as binding and unbinding of motors to or from the
microtubule. For different values of the rate of kinesin association to the
cargo, the run-lengths show prominent peaks as the crowding density is
changed. However, overall, the run-length decreases significantly due to
binding of motors to the microtubule, a process that increases the effective
crowding density. As a consequence of cargo’s ability to associate kinesin,
the association-time of the cargo to the microtubule increases with the
increase in the crowding density. The velocity of the cargo, on the other
hand, decreases with the crowding density and it remains approximately
unchanged with the kinesin-association rate of the cargo. Incorporating the
processive motion in the mathematical model would add another level of
complexity which can be a subject of future studies.

##  Appendix A Model 1

In the matrix form the differential equations (4) and (5) appear as

$\displaystyle\frac{d}{dt}{\bf G}(\gamma,t)={\bf R}{\bf G}(\gamma,t),$ (17)

where

$\displaystyle{\bf G}(\gamma,t)={\begin{pmatrix}\tilde{P}(\gamma,t)\\\
\tilde{Q}(\gamma,t)\\\ \end{pmatrix}}\ {\rm and}\ $ $\displaystyle{\bf
R}={\begin{pmatrix}&(1-r_{m})\gamma-1&r_{an}\\\ &r_{m}\gamma&-r_{an}\\\
\end{pmatrix}.}$ (18)

Here ${\bf R}$ is a transition matrix. For $\gamma=1$, the sum of all the
elements in a column is $0$. One may find out the solutions of these equations
upon finding the eigenvalues and eigenvectors of ${\bf R}$. The eigenvectors
corresponding to the eigenvalues $\lambda_{\pm}$ are, respectively,

$\displaystyle{\begin{pmatrix}&1\\\
&\frac{\lambda_{+}+1-(1-r_{m})\gamma}{r_{an}}\\\ \end{pmatrix}}\ \ {\rm and}\
\ {\begin{pmatrix}&1\\\ &\frac{\lambda_{-}+1-(1-r_{m})\gamma}{r_{an}}\\\
\end{pmatrix}},$ (19)

where $\lambda_{+,-}=\frac{1}{2}[-(r_{an}+1-(1-r_{m})\gamma)\pm A]$ with
$A=\sqrt{(r_{an}+1-(1-r_{m})\gamma)^{2}-4r_{an}(1-\gamma)}$. The solutions for
the generating functions are

$\displaystyle{\begin{pmatrix}&\tilde{P}(\gamma,t)\\\ \
&\tilde{Q}(\gamma,t)\\\
\end{pmatrix}}=c_{1}e^{\lambda_{+}t}{\begin{pmatrix}&1\\\
&\frac{\lambda_{+}+1-(1-r_{m})\gamma}{r_{an}}\\\ \end{pmatrix}}+$
$\displaystyle c_{2}e^{\lambda_{-}t}{\begin{pmatrix}&1\\\
&\frac{\lambda_{-}+1-(1-r_{m})\gamma}{r_{an}}\\\ \end{pmatrix}},$ (20)

where $c_{1},\ c_{2}$ are integration constants. We consider the initial
conditions $P(i,t=0)=Q(i,t=0)=1/2$ for $i=0$. Using these conditions, we find

$\displaystyle
c_{1}=\frac{1}{2A}\left[r_{an}-\lambda_{-}-1+(1-r_{m})\gamma\right]\ \ \ {\rm
and}$ $\displaystyle
c_{2}=\frac{1}{2}-c_{1}=\frac{1}{2A}\left[A-r_{an}+\lambda_{-}+1-(1-r_{m})\gamma\right].$
(21)

Since in the large time limit, the solutions are governed by the largest
eigenvalue, we have

$\displaystyle\tilde{P}(\gamma,t)+\tilde{Q}(\gamma,t)\approx
c_{1}e^{\lambda_{+}t}\left[1+\frac{\lambda_{+}+1-(1-r_{m})\gamma}{r_{an}}\right].$

Upon taking derivatives of $\tilde{P}(\gamma,t)+\tilde{Q}(\gamma,t)$ with
respect to $\gamma$, we have

$\displaystyle\langle
i\rangle=\left[\gamma\left(\frac{d\tilde{P}}{d\gamma}+\frac{d\tilde{Q}}{d\gamma}\right)\right]_{\gamma=1}$
$\displaystyle=\bigg{\\{}\gamma\frac{dc_{1}}{d\gamma}e^{\lambda_{+}t}\big{(}1+\frac{\lambda_{+}+1-(1-r_{m})\gamma}{r_{an}}\big{)}\bigg{\\}}_{\gamma=1}+$
$\displaystyle\bigg{\\{}\gamma
c_{1}e^{\lambda_{+}t}\frac{d\lambda_{+}}{d\gamma}t\big{(}1+\frac{\lambda_{+}+1-(1-r_{m})\gamma}{r_{an}}\big{)}\bigg{\\}}_{\gamma=1}+$
$\displaystyle\bigg{\\{}\gamma
c_{1}e^{\lambda_{+}t}\frac{1}{r_{an}}\big{(}\frac{d\lambda_{+}}{d\gamma}-(1-r_{m})\big{)}\bigg{\\}}_{\gamma=1}.$
(23)

In the large time limit, we finally have

$\displaystyle\langle i\rangle/t=\frac{d\lambda_{+}}{d\gamma}\mid_{\gamma=1}.$
(24)

Using

$\displaystyle\frac{d\lambda_{+}}{d\gamma}\bigg{|}_{\gamma=1}=\left[\frac{1-r_{m}}{2}+\frac{1}{2}\frac{dA}{d\gamma}\right]_{\gamma=1},$
(25)

where

$\displaystyle\frac{dA}{d\gamma}\bigg{|}_{\gamma=1}=\frac{r_{\rm an}+r_{\rm
an}r_{m}-r_{m}+r_{m}^{2}}{r_{\rm an}+r_{m}},$ (26)

we have

$\displaystyle
v=\frac{d\lambda_{+}}{d\gamma}\mid_{\gamma=1}=\frac{r_{an}}{r_{an}+r_{m}}.$
(27)

Figure (13) shows plots of velocity obtained mathematically and through
simulations.

Figure 13: Velocity of the cargo as a function of $r_{an}$ with $r_{m}=1/2$.

##  Appendix B Model 2

###  2.1 $m=3$

A cargo that can bind at the most three kinesins can be in four possible
states, namely, bound to one, two or three kinesins or not bound to any
kinesin. Possible configurations of two neighbouring sites can be of $P$ type
or $Q$ type depending on whether the site in front of the cargo is occupied by
a free kinesin or empty. For example, $P_{n}(i,t)$ ($n=1,\ 2,\ {\rm or}\ 3$)
indicates the probability of a configuration where a cargo, bound to $n$
number of kinesins, is present at the $i$-th site at time $t$ while the
$(i+1)$-th site is empty. Similarly, $Q_{n}(i,t)$ ($n=1,\ 2,\ {\rm or}\ 3$)
represents the probability of a configuration where a cargo, bound to $n$
number of kinesins, is present at the $i$-th site at time $t$ while the
$(i+1)$-th site is occupied by a free kinesin. The change in these
probabilities with time is described by the equations

$\displaystyle\frac{d}{dt}P_{3}(i)=(1-r_{m})P_{3}(i-1)-P_{3}(i)+r_{an}Q_{2}(i)-\omega_{d}P_{3}(i)+\
\omega_{a}P_{2}(i),$ (28)
$\displaystyle\frac{d}{dt}P_{2}(i)=(1-r_{m})P_{2}(i-1)-P_{2}(i)+r_{an}Q_{1}(i)+\omega_{d}P_{3}(i)+\omega_{a}P_{1}(i)-(\omega_{a}+\omega_{d})P_{2}(i),$
(29)
$\displaystyle\frac{d}{dt}P_{1}(i)=(1-r_{m})P_{1}(i-1)-P_{1}(i)+\omega_{d}P_{2}(i)-(\omega_{a}+\omega_{d})P_{1}(i),$
(30)
$\displaystyle\frac{d}{dt}Q_{3}(i)=r_{m}P_{3}(i-1)+\omega_{a}Q_{2}(i)-\omega_{d}Q_{3}(i),$
(31)
$\displaystyle\frac{d}{dt}Q_{2}(i)=r_{m}P_{2}(i-1)+\omega_{a}Q_{1}(i)+\omega_{d}Q_{3}(i)-(\omega_{d}+\omega_{a})Q_{2}(i)-r_{\rm
an}Q_{2}(i),\ \ {\rm and}$ (32)
$\displaystyle\frac{d}{dt}Q_{1}(i)=r_{m}P_{1}(i-1)-r_{an}Q_{1}(i)+\omega_{d}Q_{2}(i)-(\omega_{a}+\omega_{d})Q_{1}(i).$
(33)

Additionally, as in $m=2$ case, we have

$\displaystyle\frac{d}{dt}P_{0}(i)=\omega_{d}P_{1}(i)\ {\rm and}\
\frac{d}{dt}Q_{0}(i)=\omega_{d}Q_{1}(i).$ (34)

Defining the generating functions as
$\tilde{P}_{n}(\gamma,t)=\sum_{i=-\infty}^{\infty}\gamma^{i}P_{n}(i,t)$ and
$\tilde{Q}_{n}(\gamma,t)=\sum_{i=-\infty}^{\infty}\gamma^{i}Q_{n}(i,t)$, we
have

$\displaystyle\frac{d}{dt}{\bf H}(\gamma,t)={\bf S}{\bf H}(\gamma,t),$ (35)

where ${\bf H}$ is a column matrix

$\displaystyle{\bf H}(\gamma,t)={\begin{pmatrix}&\tilde{P}_{3}(\gamma,t)\\\
&\tilde{P}_{2}(\gamma,t)\\\ &\tilde{P}_{1}(\gamma,t)\\\
&\tilde{Q}_{3}(\gamma,t)\\\ &\tilde{Q}_{2}(\gamma,t)\\\
&\tilde{Q}_{1}(\gamma,t)\\\ \end{pmatrix}},$ (36)

and ${\bf S}$ is a $6\times 6$ matrix

$\displaystyle{\bf
S}={\begin{pmatrix}&(1-r_{m})\gamma-\omega_{d}-1&\omega_{a}\ &0\ &0\
&r_{an}&0\\\
&\omega_{d}&(1-r_{m})\gamma-\omega_{a}-\omega_{d}-1&\omega_{a}&0&0&r_{\rm
an}\\\ &0&\omega_{d}&(1-r_{m})\gamma-\omega_{a}-\omega_{d}-1&0&0&0\\\
&r_{m}\gamma&0&0&-\omega_{d}&\omega_{a}&0\\\
&0&r_{m}\gamma&0&\omega_{d}&-\Omega&\omega_{a}\\\
&0&0&r_{m}\gamma&0&\omega_{d}&-\Omega\\\ \end{pmatrix}},$ (37)

where $\Omega=(\omega_{a}+\omega_{d}+r_{\rm an}).$ As in case of $m=2$, here
again the variation in the run-length is governed by the quantity $t_{\rm
assoc}\frac{d\lambda_{l}}{d\gamma}\mid_{\gamma=1}$ where $\lambda_{l}$ is the
largest eigenvalue of matrix $S$. Figure (14) shows the variation in the
$t_{\rm assoc}\frac{d\lambda_{l}}{d\gamma}\mid_{\gamma=1}$ with $r_{m}$.

Figure 14:  $y$-axis represents the product of two factors that dominate the
nature of the run-length. $t_{\rm assoc}$ is the association time of the cargo
and $\frac{d\lambda_{l}}{d\gamma}$ arises while computing the average velocity
of the cargo (see the text). This plot is for $m=3$ with
$\omega_{a}=\omega_{d}=0.01$.

##  Appendix C Processive movement of free kinesins

Descriptions of parameters used in simulations are provided in table 1.

$\Delta t$ | $d/v[s]$ |  discrete time step, $v$ [nm $s^{-1}$] - velocity of free kinesin/cargo, $d$ [nm] - length of the tubulin dimer (lattice spacing)
---|---|---
$\omega_{a,kin}$ | $\bar{\omega}_{a,kin}d\Delta t$ |  dimensionless kinesin attachment rate
$\omega_{d,kin}$ | $\bar{\omega}_{d,kin}\Delta t$ |  dimensionless kinesin detachment rate
$\omega_{a}$ | $-$ |  dimensionless kinesin attachment rate to cargo
$\omega_{d}$ | $-$ |  dimensionless kinesin detachment rate from cargo
$r_{\rm an}$ | $-$ |  dimensionless association rate of free kinesin to cargo
Table 1: Description of parameters used in simulations in section 2.3. A
tubulin dimer is the basic subunit of a microtubule. $\bar{\omega}_{a,kin}$
and $\bar{\omega}_{d,kin}$. $\bar{\omega}_{a,kin}$ denotes the attachement
rate of a free kinesin to the microtubule per unit length per unit time [$\mu
m^{-1}\ s^{-1}$] and $\bar{\omega}_{d,kin}$ denotes the detachment rate of a
free kinesin from the microtubule per unit time [$s^{-1}$].

Figure (15) shows how the run-length varies with the crowding density in the
three cases - Processive movement of free kinesins and (i) binding of kinesins
to the microtubule at rate $\omega_{a,{\rm kin}}$, (ii) unbinding of kinesins
from the microtubule at rate $\omega_{d,{\rm kin}}$, and (iii) no binding or
unbinding of kinesins to or from the microtubule.

Figure 15:  $\omega_{a}=\omega_{d}=0.05$, $r_{\rm an}=0.4$, $p=0.1$, and
$\beta=0.6$. The total number of lattice sites is $2000$. Run-lengths are
obtained upon averaging over $4000$ samples. In case of only attachment,
$\omega_{a,{\rm kin}}=0.01$. In case of only detachment, $\omega_{d,{\rm
kin}}=0.1$.

## References

  * (1) J. Howard, Mechanics of Motor Proteins and the Cytoskeleton (Sinauer Associates, Massachusetts, 2001).
  * (2) M. Schliwa, Molecular Motors (Wiley-VCH, Germany, 2003).
  * (3) J. Howard, A. J. Hudspeth, and R. D. Vale, Nature 342, 154 (1989).
  * (4) K. Visscher, M. J. Schnitzer, and S. M. Block, Nature 400, 184 (1999).
  * (5) N. J. Carter, and R. A. Cross, Nature 435, 308 (2005).
  * (6) S. P. Gross, M. A. Welte, S. M. Block, and E. F. Wieschaus, J. Cell Biol. 156, 715 (2002).
  * (7) D. B. Hill, M. J. Plaza, K. Bonin, and G. Holzwarth, Eur. Biophys. J. 33, 623 (2004).
  * (8) V. Levi, A. S. Serpinskaya, E. Gratton, and V. Gelfand, Biophys. J. 90, 318 (2006).
  * (9) K. J. Böhm, R. Stracke, P. Mühlig, and E. Unger, Nanotechnology 12, 238 (2001).
  * (10) S. Klumpp, and R. Lipowsky, Proc. Natl. Acad. Sci. U. S. A. 102, 17284 (2005).
  * (11) J. Beeg, S. Klumpp, R. Dimova, R. S. Graciá, E. Unger, and R. Lipowsky, Biophys. J. 94, 532 (2008).
  * (12) M. Vershinin, B. C. Carter, D. S. Razafsky, and S. P. Gross, Proc. Natl. Acad. Sci. U. S. A. 104, 87 (2007).
  * (13) M. E. Schneider, et al, J. Neurosci. 26, 10243 (2006).
  * (14) P. I. Zhuravlev, B. S. Der, and G. A. Papoian, Biophys. J. 98, 1439 (2010).
  * (15) C. Leduc et al., Proc. Natl. Acad. Sci. U. S. A. 109, 6100 (2012).
  * (16) V. Bormuth, et. al. Biophys. J. 103, L4-L6 (2012); M. Bugiel, E. Böhl and E. Schäffer, Biophys. J. 108, 2019 (2015).
  * (17) C. T. MacDonald, J. H. Gibbs, and A. C. Pipkin, Biopolymers 6, 1 (1968).
  * (18) M. R. Evans, D. P. Foster, C. Godrèche, and D. Mukamel, Phys. Rev. Lett. 74, 208 (1995).
  * (19) A. Parmeggiani, T. Franosch, and E. Frey, Phys. Rev. Lett. 90, 086601 (2003).
  * (20) Y. Chai, S. Klumpp, M. J. I. Müller, and R. Lipowsky, Phys. Rev. E 80, 041928 (2009).
  * (21) A. I. Curatolo, M. R. Evans, Y. Kafri, and J. Tailleur, J. Phys. A 49, 095601 (2016).
  * (22) I. Pinkoviezky and N. S. Gov, Phys. Rev. Letts. 118, 018102 (2017).
  * (23) P. Wilke, E. Reithmann, and E. Frey Phys. Rev. X 8, 031063 (2018).
  * (24) A. Seitz, and T. Surrey, EMBO J. 25, 267 (2006).
  * (25) L. Conway, et. al., Proc. Natl. Acad. Sci. U. S. A. 109, 20814 (2012).
  * (26) L. Conway, and J. Ross, Comm. and Integ. Biol. 6, e25387.
  * (27) Sutapa Mukherji, Phys. Rev. E 77, 051916 (2008).


figurem tablem

Direct validation of dune instability theory
— Supplementary Information —

Ping Lü1,∗, Clément Narteau2,∗, Zhibao Dong1,

Philippe Claudin3, Sébastien Rodriguez2, Zhishan An4,

Laura Fernandez-Cascales2, Cyril Gadal2, Sylvain Courrech du Pont5

  * 1

School of Geography and Tourism, Shaanxi Normal University, 620 Chang’an West
Avenue, Xi’an, Shaanxi 710119, China.

  * 2

Université de Paris, Institut de physique du globe de Paris, CNRS, F-75005
Paris, France.

  * 3

Physique et Mécanique des Milieux Hétérogènes, UMR 7636 CNRS, ESPCI PSL
Research Univ, Sorbonne Univ, Université de Paris, 10 rue Vauquelin, 75005
Paris, France.

  * 4

Northwest Institute of Eco-Environment and Resources Donggang West Road 320,
Lanzhou, Gansu Province 730000, China.

  * 5

Laboratoire Matière et Système Complexes, Université de Paris. UMR 7057 CNRS,
Bâtiment Condorcet, 10 rue Alice Domon et Léonie Duquet, 75205 Paris Cedex 13,
France.

∗To whom correspondence should be addressed. E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>

###### Contents

  1. 1 Supplementary Note 1  Wind data
  2. 2 Supplementary Note 2  Grain size distribution
  3. 3 Supplementary Note 3 Transport properties on a flat sand bed
  4. 4 Supplementary Note 4 Estimating the saturation length $l_{\rm sat}$ in the field
  5. 5 Supplementary Note 5 Upwind velocity shift at dune crests and troughs
    1. 5.1 The inner and outer layers
    2. 5.2 Measuring the upwind velocity shift within the inner layer
    3. 5.3 From the upwind velocity shift to the hydrodynamic parameters $A$ and $B$
  6. 6 Supplementary Note 6 Dispersion diagram from successive topographic surveys
    1. 6.1 Aeolian dune formation as a linear instability
    2. 6.2 The transition from the linear to the non-linear phases of dune growth
    3. 6.3 A transport time scale for dune growth under variable wind strength
    4. 6.4 Topographic surveys during experimental dune growth
    5. 6.5 Spectral analysis, amplitude and phase
    6. 6.6 Dispersion diagram of the growth rate
    7. 6.7 Dispersion diagram of the phase velocity

##  1 Supplementary Note 1
Wind data

We installed a 2 m high wind tower in the center of our new experimental site
and collected the wind data of the local airport located 10 km east. We use
the local tower for the measurements of the transport threshold and the
saturation length (see Figs. 2 and 3 of the main manuscript). We check the
consistency between the local wind data and those collected at the airport.
For the long-term experiment, we only use the wind data of the local airport.

Figure S1: Comparison of wind speed data. Wind shear velocity derived from
temporary wind speed measurements made during the threshold experiment (blue)
and simultaneous measurements obtained on the local (red) and the airport wind
towers (yellow). There is a general agreement between all these data sets.

Fig. S1 shows the shear velocity derived from a local wind measurement, the 2
m high wind tower and the airport meteorological tower during the transport
threshold experiment (see Fig. 2 of the main manuscript). All these data are
consistent with each other. As a consequence, the threshold shear stress
derived from the local measurement can be extrapolated to other wind data to
compute sand fluxes (Sec. 3).

##  2 Supplementary Note 2
Grain size distribution

The landscape-scale experiment site is located close to the oasis city of
Shapotu at 8 km from the Yellow River in the Tengger Desert, which covers an
area of about 36,700 km2 in the northwest part of the Zhongwei County in the
Ningxia Hui Autonomous Region of the People’s Republic of China (37° 31´ N,
105° E). This desert is characterized by a lognormal grain size distribution
with a mean value $d$=190$\leavevmode\nobreak\ \mu$m (Fig. S2).

Figure S2: Grain size distribution in the Tengger Desert. The line is the best
fit using a lognormal distribution with a mean value $d=190\;\mu{\rm m}\approx
10^{-3.72}\;{\rm m}$.

Variable | Units | Value
---|---|---
Acceleration of gravity $g$ | ${\rm m\,s}^{-2}$ | 9.81
Grain size $d$ | m | 190$\times 10^{-6}$
Air density $\rho_{\rm f}$ | ${\rm kg\,m}^{-3}$ | 1.29
Grain density $\rho_{\rm s}$ | ${\rm kg\,m}^{-3}$ | 2.55$\times 10^{3}$
Aerodynamic roughness $z_{\rm s}$ | m | 10-3
von-Kármán constant $\kappa$ | $\emptyset$ | 0.4
Shear velocity and sand flux on a flat sand bed
Threshold shear velocity $u_{\rm th}\;\;\;\;\;\;\;\;\;\;\;$ | m s-1 | 0.19
Mean shear velocity $\langle u_{*}\rangle$ | m s-1 | 0.29
$\langle u_{*}\rangle/u_{\rm th}$ | $\emptyset$ | 1.5
${\rm DP}=\langle\|\overrightarrow{Q}\|\rangle$ | m${}^{2}\;$yr-1 | 18.4
${\rm RDP}=\|\langle\overrightarrow{Q}\rangle\|$ | m${}^{2}\;$yr-1 |  05.7
RDP/DP | $\emptyset$ | 0.32
${\rm RDD}$ | mod 360° | 266.6

Table S1:  Shear velocity and sand fluxes derived from the airport wind data
from January 1, 2013 to October 31, 2017. The wind and sand flux roses are
shown in Fig. 1B of the main manuscript. See text and Eqs. 1-9 for the
description of all the variables. The resultant drift direction (${\rm RDD}$)
is measured counterclockwise from East.

##  3 Supplementary Note 3
Transport properties on a flat sand bed

Wind data are used to predict sand flux properties on a flat sand bed. Tab. S1
shows the results obtained using the formalism that follows.

Wind measurements provide the wind speed $u_{i}$ and direction $\vec{x}_{i}$
at different times $t_{i},\;i\in[1;N]$. For each time step $i$, the shear
velocity writes

$u_{*}^{i}=\displaystyle\frac{u_{i}\kappa}{\log(z/z_{s})},$ (1)

where $z$ is the height at which the wind velocity $u_{i}$ has been measured
and $\kappa$ the von-Kármán constant. Instead of the geometric roughness that
depends only on grain size, we consider here the aerodynamic roughness $z_{\rm
s}$ that accounts for the height of the transport layer in which saltating
grains modify the vertical wind velocity profile. The value of the threshold
shear velocity for motion inception is determined using the formula calibrated
by Iversen and Rasmussen[1]

$u_{\rm th}=0.1\displaystyle\sqrt{\frac{\rho_{\rm s}}{\rho_{\rm f}}gd}.$ (2)

Using the gravitational acceleration $g$, the grain to fluid density ratio
$\rho_{\rm s}/\rho_{\rm f}\simeq 2.05\times 10^{3}$ and the grain diameter
$d=190\,\mu{\rm m}$, we find $\rm u_{\rm th}=0.19\,m\,s^{-1}$, which
corresponds to a threshold wind speed of $u_{10}=4.4\,{\rm m}\,{\rm s}^{-1}$
ten meters above the ground. It is close to the values we measured in the
field, which are $u_{\rm th}=0.23\pm 0.0.4$ and $u_{10}=5.3\pm 0.92\,{\rm
m}\,{\rm s}^{-1}$.

For each time step $i$, the saturated sand flux $\overrightarrow{Q_{i}}$ on a
flat sand bed is computed from the relationship proposed by Ungar and Haff[2]
and calibrated by Durán et al.[3]

$Q_{\rm sat}(u_{*})=\begin{cases}25\,\displaystyle\frac{\rho_{\rm
f}}{\rho_{\rm s}}\sqrt{\displaystyle\frac{d}{g}}\left(u_{*}^{2}-u_{\rm
th}^{2}\right)&\text{{\small for} $u_{*}>u_{\rm th}$},\\\ 0&\text{{\small
else}}.\end{cases}$ (3)

In this formula, the prefactor takes into account a dune compactness of $0.6$.

From the individual saturated sand flux vectors $\overrightarrow{Q_{i}}$, we
estimate the mean sand flux vector on a flat erodible bed

$\langle\overrightarrow{Q}\rangle=\displaystyle{\displaystyle\sum_{i=2}^{N}\overrightarrow{Q_{i}}\delta
t_{i}}\bigg{/}{\displaystyle\sum_{i=2}^{N}\delta t_{i}},$ (4)

where $\delta t_{i}=t_{i}-t_{i-1}.$ The norm of the mean sand flux is usually
called the resultant drift potential:

${\rm RDP}=\|\langle\overrightarrow{Q}\rangle\|.$ (5)

This quantity is highly dependent on the wind regime. Since it is a vectorial
sum, the contributions of winds from opposite directions cancel each other
out. For the entire time period, we also calculate the drift potential,

${\rm
DP}=\displaystyle{\displaystyle\sum_{i=2}^{N}\left\|\overrightarrow{Q_{i}}\right\|\delta
t_{i}}\bigg{/}{\displaystyle\sum_{i=2}^{N}\delta t_{i}}.$ (6)

Unlike the resultant drift potential, this mean sand flux does not take into
account the orientation of the individual sand fluxes computed from the
successive wind measurements[4].

The ratio $\rm RDP/DP$ is a non-dimensional parameter, which is often used to
characterize the directional variability of the wind regimes[5, 6]: $\rm
RDP/DP$ $\to 1$ indicates that sediment transport tends to be unidirectional;
$\rm RDP/DP$ $\to 0$ indicates that most of the transport components cancel
each other. Finally, RDD is the resultant drift direction, i.e., the direction
of $\langle\overrightarrow{Q}\rangle$.

The mean shear velocity $\langle u_{*}\rangle$ is defined as the shear
velocity averaged over the transport periods. i.e. when $Q_{\rm sat}>0$. Using
the Heaviside function $H_{u_{\rm th}}$ defined as

$H_{u_{\rm th}}=\begin{cases}1&\text{{\small for} $u_{*}>u_{\rm th}$},\\\
0&\text{{\small else}}.\end{cases}$ (7)

the mean shear velocity can be defined directly from the shear velocity

$\langle u_{*}\rangle={\displaystyle\sum_{i=2}^{N}H_{u_{\rm
th}}^{i}u_{*}^{i}}\bigg{/}{\displaystyle\sum_{i=2}^{N}H_{u_{\rm th}}^{i}},$
(8)

or from the integrated flux using the inverse function $Q_{\rm sat}^{-1}$ of
the transport law (Eq. 3)

$\langle u_{*}\rangle=Q_{\rm sat}^{-1}\left({\rm
DP}\times{\displaystyle\sum_{i=2}^{N}\delta
t_{i}}\bigg{/}{\displaystyle\sum_{i=2}^{N}H_{u_{\rm th}}^{i}\delta
t_{i}}\right).$ (9)

These two estimations of $\langle u_{*}\rangle$ are close to each other
considering wind data from the Tengger Desert.

##  4 Supplementary Note 4
Estimating the saturation length $l_{\rm sat}$ in the field

Let us consider an infinite flat granular bed under a unidirectional wind in a
statistically steady state. Eventually, the transport rate reaches an
equilibrium state due to the negative feedback between the density of moving
grains and the strength of the flow, which determines the saturated flux
$Q_{\rm sat}$ (Eq. 3). We consider now a situation for which the flow and the
sediment flux is non-homogeneous or unsteady in space or time. The actual flux
$q$ does not immediately adjust to the local value of the shear stress[7, 8,
9]. It needs some space and time to reach its equilibrium $Q_{\rm sat}$-value
(Fig. S3). Over bedforms, the transport is never far from its saturated state,
so it can be expressed by a first-order linear relaxation in both space and
time

$t_{\rm sat}\,\dfrac{\partial q}{\partial t}+l_{\rm sat}\,\dfrac{\partial
q}{\partial x}=q_{\rm sat}-q,$ (10)

where $l_{\rm sat}$ and $t_{\rm sat}$ are the saturation length and the
saturation time, respectively[8, 10, 11]. The $t_{\rm sat}$-value is usually
much smaller ($\approx 1$ s) than the characteristic time scale for the
evolution of the bed ($\approx 10^{5}$ s). Such a separation in scale
justifies the simplifying assumption that the fluid flow and sediment
transport can be considered and computed as if the bed was fixed[3].
Neglecting $t_{\rm sat}$, the saturation transient described by Eq. 10 has
been successfully applied to the description of dune formation[9, 11, 10, 12,
13]. Under the configuration shown in Fig. S3, we have

$q(x)=Q_{\rm sat}\left(1-\exp\left(\dfrac{x-x_{0}}{l_{\rm
sat}}\right)\right).$ (11)

In practice, the saturation length $l_{\rm sat}$ scales as $(\rho_{\rm
s}/\rho_{\rm f})d$ the distance needed for one grain to be accelerated up to
the wind velocity[10, 14, 15, 16]:

$l_{\rm sat}\approx 2.2\,\dfrac{\rho_{\rm s}}{\rho_{\rm f}}d,$ (12)

Figure S3:  A schematic representation of the saturation length experiment.
The saturation length may be described as the relaxation length of the sand
flux $q$ toward its saturated value $Q_{\rm sat}$. The experiment set-up is a
flat sand bed preceded by a flat non-erodible bed. When the wind blowing over
this rigid surface reaches the erodible bed, the sand flux relaxes
exponentially towards its saturation value $Q_{\rm sat}$ over a characteristic
distance $l_{\rm sat}$ in the direction of the flow.

As shown in Fig. S4 and in Fig. 3 of the main manuscript, we determine the
value of $l_{\rm sat}$ from the evolution of the topography of a rectangular
sand pile placed downstream of a non-erodible bed composed of gravels. We thus
have the experimental setup that best reflects the theoretical conditions
presented in Fig. S3. The rectangular sand pile has a length of 12 m and a
width of 3 m with its main axis aligned in the northwest-southeast direction,
parallel to the orientation of the prevailing wind (Fig. S4a). The initial
sand pile was prepared and scanned on April 22, 2014 (Fig. S4b). A storm
occurred on April 24 with winds from the north-northwest and irregular wind
speeds reaching 15 m s-1 at a height of 2 m. The sand pile was scanned again
on April 25 (Fig. S4c).

Considering the ideal assumption of our initially flat surface under a steady
wind of constant direction, strength and transport rate $q$, Eq. 11 can be
combined with the equation of conservation of mass to estimate the erosion
rate

$\dfrac{\partial h(x)}{\partial t}=-\dfrac{\partial q(x)}{\partial
x}=-\dfrac{Q_{\rm sat}}{l_{\rm sat}}\exp\left(\dfrac{x-x_{0}}{l_{\rm
sat}}\right).$ (13)

This equation can be integrated to estimate the net erosion over a given time
period. In this case the exponential regime is expected to hold. However, it
cannot be observed from the start of the erodible bed for several reasons,
both theoretical and experimental. The main reasons are related to (1) the
natural variability of wind speed and direction and (2) the development of a
discontinuity in the topographic profile between the non-erodible and the
erodible beds. In addition, Eq. 13 neither takes into account the possible
dependence of $l_{\rm sat}$ on wind speed[17], nor the spatial shift
associated with the establishment of a transport layer when the sand flux
starts from zero[16].

Using the field data and despite the number of simplifying assumptions, we
study the difference in height $\Delta H$ between the two surface elevations
to look for zones where

$\Delta H\sim-\exp(x/l_{\rm sat}).$ (14)

In practice, we use 2D elevation profiles aligned with the primary sand
transport direction during the storm. It forms an angle of 22° with the
orientation of the initial sand pile. This angle is estimated from the sand
flux rose as well as from the orientation of ripples and accretion mounds
(Fig. S4c).

The slope of the exponential decay in Eq. 14 gives the value of $l_{\rm sat}$
(Fig. 3D of the main manuscript). We find a value of 0.95 m from our field
data, a value that can be directly compared to the 0.83 m predicted by Eq. 12
and the parameters given in Tab. S1.

Figure S4:  Measuring the saturation length $l_{\rm sat}$ in the field. (a)
Setting up the $l_{\rm sat}$ experiment between the two sand piles of the
elongating dune experiment in a triangular area of bare sand within a flat bed
armored with gravels. In this area, a rectangular sand pile 12 m long and 3 m
wide is placed downstream of the non-erodible gravel bed. The long axis is
aligned along the northwest-southeast direction, the orientation of the
prevailing wind. The sand bed is flattened so that there is no vertical step
between the erodible and non-erodible zones. (b) Elevation of the rectangular
sand pile on April 22, 2014 before aeolian transport. (c) Elevation of the
rectangular sand pile on April 25, 2014 after the passage of a northwestern
depression (see flux rose in Fig. 3C of the main manuscript). $H$ is the
absolute elevation from a base level. $H^{\prime}$ is the difference of
elevation with a smoothed surface obtained using a sliding window with a
radius of 20 cm. The inset in c show the transect along which the elevation
profile used to estimate $l_{\rm sat}$ has been extracted (see Fig. 3B of the
main manuscript). It is parallel to the sand deposits that form downwind of
topographic obstacles and perpendicular to aeolian ripples.

##  5 Supplementary Note 5
Upwind velocity shift at dune crests and troughs

Figure S5:  Experimental set-up to estimate the upwind velocity shift at the
crests and troughs of incipient dunes. (a) The mobile anemometer mast, with
anemometers located at heights of 4 and 12 cm in the inner layer and at
heights of 50 and 100 cm in the outer layer. Arrows show the mean wind
direction and the measurement path. (b) The anemometer mast and the reference
anemometer. (c) A measurement path perpendicular to the crest orientation of
the bed instability. (d) Setting a new measurement for 10 min with a sampling
frequency of 1 Hz.

###  5.1 The inner and outer layers

Flows that are topographically forced by obstacles (e.g., hills or sand dunes)
accelerates on the upwind slopes and decelerates on the downwind slopes. Then,
in order to study sediment transport over a dune, we should first describe
some properties of the turbulent flow over an undulating topography.
Conceptually, as proposed by Jackson and Hunt (1975)[18] in the limit of small
amplitudes bedforms, the turbulent flow over a sinuous bed elevation profile
of wavelength $\lambda$ can be decomposed into two layers:

  * •

The outer layer is the external region (supposed infinite) where the pressure
gradient set up by the topography is balanced by the inertial forces. The
streamlines follow the topography. At a given height, the wind speed is
maximum (minimum) above the top (bottom) of the topography. The amplitude of
the perturbation vanishes on a characteristic height that varies according to
the wavelength of the topography. Then, sufficiently far above the obstacle,
the wind speed is finally equal to the undisturbed wind speed in the absence
of topography (see Eq. 1).

  * •

The inner layer is a zone in which the longitudinal pressure gradient exerted
by the fluid is compensated by the Reynolds shear stress induced by the
turbulent motions at the surface of the bed. Then, the pressure gradient is in
phase quadrature with the topography and is maximum where the stoss slope is
steepest. Hence, there is an upwind velocity shift within the inner layer. The
characteristic thickness $l$ of the inner layer depends on both the wavelength
of the bed elevation profile and on the aerodynamic roughness $z_{\rm s}$ so
that[18, 19]

$\dfrac{l}{\lambda}\log^{2}\left(\dfrac{l}{z_{s}}\right)=2\kappa^{2}.$ (15)

For typical value of the aerodynamic roughness ($z_{\rm s}<10^{-2}$ m) and
wavelengths of hundreds of meters, the inner layer is always confined in a
meter scale envelop above the topography. For $\lambda=20\;{\rm m}$ and
$z_{\rm s}=10^{-3}\;{\rm m}$, common values observed during the development of
aeolian bedforms under high wind speed, the inner layer is less than 20 cm.

###  5.2 Measuring the upwind velocity shift within the inner layer

A sufficiently strong wind and a stable wind orientation are essential for the
quality of the measurements. The sequence of measurements consists in moving
an anemoter mast upwind along a given profile of known elevation. The entire
procedure to estimate the upwind velocity shift on dunes is described in full
details in Claudin et al. (2013)[20]. We follow here the same procedure with
the specifications detailed below.

We performed our measurements on low dunes with sinusoidal shape. The density
of measurements is high near the crests and troughs ($\approx 20\;{\rm cm}$)
but lower in the steeper sections ($\approx 1\;{\rm m}$). Each measurement
lasts 10 min with a sampling frequency of 1 Hz. The mobile anemometer mast has
sonic anemometers located at heights of 4 and 12 cm and cup anemometers at
heights of 50 and 100 cm (Fig. S5a). Given the natural variability in wind
speed between two measurements, the wind speeds measured along the mobile mast
are normalized by the wind speed measured at one meter high by a reference cup
anemometer located at the top of a larger dune in the vicinity (Fig. S5b). The
normalized wind speeds along the profile can then reveal how sensitive is the
flow to topography at different heights. As shown in Fig. 4B of the main
manuscript, in the limit of small sinusoidal oscillations, the wind profiles
at all heights reflect the topography, i.e., the spatial variation of the wind
speed exhibits the same wavelength as the bedforms. Most importantly for our
present purposes, the flow is in phase with the topography at heights of 50
and 100 cm and in phase advance at heights of 4 and 12 cm. It indicates that
the inner layer has a thickness between 12 and 50 cm, a range of value that
covers entirely the values given by Eq. 15 for aeolian dune systems on Earth.

###  5.3 From the upwind velocity shift to the hydrodynamic parameters $A$ and
$B$

On several occasions during our field campaigns, we have measured wind
velocity in the inner and outer layers on a succession of crests and troughs
to estimate the phase shift over more than a wavelength of the dune pattern
(Fig. S5c,d). Here, we focus only on individual dune crests or troughs, which
are approximated by

$H(x)=H_{\rm ref}+h\cos(kx),$ (16)

where $k$ is the wave number (wavelength $\lambda=2\pi/k$), $h$ the amplitude.
The position of the dune crests and troughs are set at $x=0$ and
$x=\lambda/2$, respectively. The arbitrary reference level $H_{\rm ref}$ is
here chosen such that $H>0$. Fitting Eq. 16 to the dune elevation data gives
the values of $k$ and $h$. The dune aspect ratio is $R=2h/\lambda=hk/\pi$.
When the aspect ratio is below $\simeq 0.05$, we expect that the (low)
perturbation calculation of the aerodynamics to be valid. A more refined
analysis would involve the computation of the Fourier transform of the dune
profile, in order to account for a whole range of wave numbers $k$.

Alike topography, the wind profile at a given height along the dune can also
be fitted by a sinusoidal function of the same wave number $k$ as for the dune
elevation:

$u_{\rm b}(x)=u_{\rm b}^{0}+\delta u_{\rm b}\cos(kx+\varphi_{\rm b}).$ (17)

As said in Sec. 5.1, the topography $Z(x)$ and the wind velocity $u_{\rm
b}(x)$ are not in phase for anemometers at heights of 4 and 12 cm: the
velocity reaches its maximum upstream of the crest. Here the phase difference
is $\varphi_{\rm b}$ so that the upwind velocity shift in the inner layer is
$\varphi_{\rm b}/k$.The two other fitting parameters are $u_{\rm b}^{0}$ and
$\delta u_{\rm b}$.

Because the logarithmic law of the wall (Eq. 1) locally holds in the inner
layer at each position $x$, the velocity can be used as a proxy to calculate
the basal shear stress with $\tau_{\rm b}\propto\rho_{\rm f}u_{\rm b}^{2}$. By
expansion to the first order, we can then write:

$\tau_{\rm b}(x)\propto\rho_{\rm f}(u_{\rm
b}^{0})^{2}\left(1+kh\left(A\cos(kx)-B\sin(kx)\right)\right),$ (18)

where $A$ and $B$ are given by

$A=2\times\frac{\delta u_{\rm b}}{u_{\rm b}^{0}}\times\dfrac{\cos(\varphi_{\rm
b})}{kh}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\mbox{and}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;B=2\times\frac{\delta
u_{\rm b}}{u_{\rm b}^{0}}\times\dfrac{\sin(\varphi_{\rm b})}{kh}\,.$ (19)

Eqs. 18 and 19 are specific to the quadratic relation between $\tau_{\rm b}$
to $u_{\rm b}$. Since sediment transport is controlled by the basal shear
stress, these two parameters $A$ and $B$ are all that are needed as
aerodynamical inputs in bedform evolution models. As described in the
introduction of the main manuscript, they are of fundamental importance for
the understanding of the mechanisms of dune growth and dune size-election (see
also Sec. 6 and Eqs. 23-24). Tab. S2 shows the results obtained in April and
November 2015 in our experiment. We observe no significant trend in the
variation of $A$ and $B$ during dune growth.

Variable |  | Units |  |  Description
---|---|---|---|---
Topography
$\lambda$ |  | m-1 |  |  dune wavelength
$k$ |  | m-1 |  |  Dune wave number
$h$ |  | m-1 |  |  Amplitude of the dune
$R$ |  | $\emptyset$ |  |  Dune aspect ratio
Flow
$\delta u_{\rm b}$ |  | ${\rm m\;s}^{-1}$ |  |  Amplitude of wind speed variation in the inner layer.
$u_{\rm b}^{0}$ |  | ${\rm m\;s}^{-1}$ |  |  Mean wind speed in the inner layer.
$\varphi_{\rm b}$ |  | o |  |  Phase shift between the flow and the topography in the inner layer.
Hydrodynamic parameters
$A$ |  | $\emptyset$ |  |  In-phase hydrodynamic parameter
$B$ |  | $\emptyset$ |  |  In-quadrature hydrodynamic parameter

| 14/04 | 15/04 | 15/04∗ | 15/04 | 16/04 | 18/04 | 03/11 | 03/11∗ | 03/11 | 15/11
---|---|---|---|---|---|---|---|---|---|---
Topography
$\lambda$ | 14.84 | 18.57 | 10.50 | 16.45 | 21.92 | 14.51 | 25.05 | 14.49 | 21.06 | 14.33
$k$ | 0.423 | 0.338 | 0.598 | 0.382 | 0.287 | 0.433 | 0.251 | 0.437 | 0.298 | 0.439
$h$ | 0.297 | 0.262 | 0.185 | 0.414 | 0.358 | 0.338 | 0.662 | 0.323 | 0.264 | 0.361
$R$ | 0.020 | 0.014 | 0.018 | 0.025 | 0.016 | 0.023 | 0.026 | 0.022 | 0.012 | 0.025
Flow
$\delta u_{\rm b}$ | 0.126 | 0.111 | 0.108 | 0.198 | 0.0624 | 0.143 | 0.306 | 0.134 | 0.132 | 0.174
$u_{\rm b}^{0}$ | 0.683 | 0.650 | 0.656 | 0.678 | 0.743 | 0.701 | 0.568 | 0.622 | 0.739 | 0.661
$\varphi_{\rm b}$ | 25.51 | 11.89 | 15.95 | 20.79 | 45.95 | 10.92 | 11.36 | 19.13 | 17.42 | 9.80
Hydrodynamic parameters
$A$ | 2.648 | 3.770 | 2.860 | 3.453 | 1.138 | 2.737 | 6.362 | 2.906 | 4.328 | 3.278
$B$ | 1.264 | 0.794 | 0.817 | 1.311 | 1.176 | 0.528 | 1.278 | 1.008 | 1.358 | 0.566

Table S2: Topography and wind measurements carried out in April and November
2015. Measurements on dune troughs are shown with a ∗ on April 4 and November
3, 2015. They are located between two dune crests along the same transect to
estimate the phase shift over more than one wavelength of the dune pattern.

##  6 Supplementary Note 6
Dispersion diagram from successive topographic surveys

Dispersion relations are plotted as growth rate $\sigma$ and phase velocity
$c$ with respect to wave number $k$. They are used in linear stability
analysis to identify the most unstable wavelength, which is likely to be
observed. We describe here how we measure the growth rate of dunes as a
function of their wave numbers $k$ from successive topographic surveys
throughout the duration of our experiment. Before, we present the rational for
studying dune formation as a linear instability in a landscape scale
experiment.

###  6.1 Aeolian dune formation as a linear instability

In what follows, the term instability characterizes the growth of a small
perturbation in a system, which is often considered to be as homogeneous and
simple as possible. Conversely, stability refers to the ability of this system
to return to its original state when perturbed. The main objective of
stability analysis is to identify the ranges of perturbation wavelengths over
which the system exhibits stable or unstable behavior. Quantitative approaches
consists of estimating the initial growth rate of each mode (wavelength)
considering that all these modes are independent of each other. Obviously, the
highest and zero growth rates are particularly important. The highest ones are
associated with the most unstable modes, which will have the greatest impact
on the system. Those at zero are neutral modes that often mark the transition
between stable and unstable regimes.

The analysis of the time and length scales of instabilities by means of
linearized equations is a standard approach in hydrodynamics and many other
branches of physics. By identifying stabilizing and destabilizing mechanisms,
these linear stability analysis reveal how they together govern the the
evolution of the system over the entire range of possible wavelengths. During
the first stage of the instability (the linear regime), the initial wavelength
$\lambda_{\rm g}$ of the perturbation is assumed to stay constant, whereas its
amplitude $h$ grows exponentially with time

$h(k_{\rm g},t)\propto\exp(\sigma(k_{\rm g})t).$ (20)

In this expression, $k_{\rm g}=2\pi/\lambda_{\rm g}$ is the wave number and
$\sigma(k_{\rm g})$ the growth rate with units of frequency. Positive and
negative growth rates correspond to unstable and stable modes, respectively.
The dispersion relation $\sigma(k)$ gives the growth rate value of the
perturbation as a function of the wave number $k$. The largest positive
$\sigma$-value corresponds to the most unstable mode $\\{k_{\rm
max},\,\lambda_{\rm max}\\}$. Zero values corresponds to neutral modes
$\\{k_{0},\,\lambda_{0}\\}$. In diffusive systems, dispersion relations are
often characterized by a transition from a stable to an unstable regime. for
an increasing wavelength. In other words, $\sigma(k)>0$ for $k<k_{0}$ and
$\sigma(k)<0$ for $k>k_{0}$. Beyond the linear stage of the instability, when
the amplitude of the initial perturbation is too high, the dispersion relation
no longer applies because the different modes interact with one another. It is
described as the non-linear stage of the instability.

Here, we study the formation of aeolian dunes as a linear instability. More
exactly, we focus on the dependence and feedback between bed forms, wind flow
and sand transport properties. The governing equations lead to the theoretical
expression of a dispersion relation as a function of the different physical
parameters of the system (see Eq. 1 of the main manuscript). Thus, the minimum
size for dunes is associated with a neutral mode and a transition from
unstable to stable regime for decreasing wavelength. The emergence and growth
of periodic dune patterns are associated with a most unstable mode which is
going to prevail within the whole dune field. Based on observation of mature
aeolian dunes in nature, previous studies have measured the wavelength and the
migration rate of dunes in order to derive values of $\lambda_{0}$,
$\lambda_{\rm max}$ and $\sigma_{\rm max}$ under various conditions.
Dispersion relations have been given less attention or even disregarded,
certainly because of the length and time scales involved in the mechanism of
aeolian dune growth. Then, a direct validation of the dune instability theory
is to investigate whether or not dispersion diagrams can be derived from field
data. Another solution is to verify whether the theoretical formalism used
with the values of the underlying physical parameters measured independently
in the field is actually capable of accurately predicting the observations.

Here, we apply this methodology to the formation of dunes in a landscape scale
experiment. Our initial system is a sand bed after flattening by a bulldozer.
In such an experimental set up, but also in all natural dune environments,
there are always heterogeneities and defects at all length scales and our
initial condition already contains all the perturbation wavelengths.

###  6.2 The transition from the linear to the non-linear phases of dune
growth

The continuous transition from the linear to the non-linear phases of dune
growth is controlled by dune aspect-ratio, which is the main control parameter
for aerodynamic non-linearities. Figs. S6a-b show the elevation and slope maps
during incipient dune growth in our experimental field from April 10, 2014 to
July 6, 2015. There is a significant change in slope maps between October 30
and November 12, 2014. Over this time interval, the local slopes not only
become steeper, they also become spatially organized, just like dune crest, to
form more regular transverse structures across the flattened area. This occurs
for a mean slope of about 0.03 (Fig. S6c) at the same time as incipient slip
faces a few centimeter high emerge (Fig. S7). The mean slope
$\langle\|\vec{\nabla}H\|\rangle$ variations can be compared to the amplitude
of the bedforms in Figs. S6b-c. Based on these observations, we set the
transition between October 30 and November 12, 2014. This transition is not
spontaneous but it is surely completed in April 2015, when the mean slope
reaches a value of about 0.07.

Figure S6:  The transition from the linear to the non-linear phases of dune
growth. (a) Elevation and (b) slope maps during incipient dune growth from
April 10, 2014 to July 6, 2015. The colormap in Fig. S6b is saturated to
distinguish the slopes of bedforms in fall 2014, around the transition from
the linear to the non-linear phases of dune growth. In non-saturated areas,
ripples can be observed. A cable crosses the experimental field from 28 April
2015. (c) The mean amplitude $2\sqrt{2}(\langle h^{2}\rangle-\langle
h\rangle^{2})^{1/2}$ of bedforms and (d) the mean slope
$\langle\|\vec{\nabla}H\|\rangle$ with respect to time. The higher mean slope
value on April 10, 2014, after the flattening, is due to bulldozer tracks that
can be seen in Fig.S6b. Colors and dashed lines are used to separate the
linear (orange) and the non-linear (blue) phases of dune growth. For
comparison with slope data, Fig. S6d is the same as Fig. 5C of the main
manuscript.

Figure S7:  Bedforms during the transition from the linear to the non-linear
phases of dune growth. (a) Incipient slip faces a few centimetres high on
November 6, 2014. (b) Periodic dune patterns on November 10, 2014. This
pictures can be compared to elevation and slope maps of October 30, 2014 and
November 12, 2014 shown in Fig. S6a and S6b.

###  6.3 A transport time scale for dune growth under variable wind strength

In order to quantify dune growth under variable wind strength, it is necessary
to define a new time scale that accounts for the intensity of transport. For
example, periods during which there is no transport should not be used as time
increments, while storm periods should contribute more to the total time.

In practice, the new time scale is defined sequentially from the wind data
using the saturated sand flux (Eq. 3) and the saturation length $l_{\rm sat}$.
Starting at $t=0$ at the flattening time, the new times write

$t_{\rm
a}(t)=\dfrac{\displaystyle\sum_{i=1}^{t}\|\overrightarrow{Q_{i}}\|\delta
t_{i}}{l_{\rm sat}^{2}}.$ (21)

This time scale is dimensionless allowing for comparison across different time
series and various time periods. Using the saturated sand flux derived from
the local wind data (Eq. 3), Fig. S8 shows the dimensionless times with
respect to time from April 10, 2014 to October 31, 2017 using the saturation
length $l_{\rm sat}=0.95\leavevmode\nobreak\ {\rm m}$ measured in the field
(Fig. 3 of the main manuscript). By definition, periods of stronger winds are
associated with steeper slopes, and vice versa. The choice of the
characteristic length for the computation of the dimensionless transport time
scale, (here $l_{\rm sat}$) is of minor significance since we consider it to
be constant.

Figure S8:  A dimensionless transport time scale for dune growth under
variable wind strength. Using the saturated sand flux derived from wind data
(Eq. 3), the red solid line shows the dimensionless transport time scale
$t_{a}$ as a function of time from April 10, 2014, when the dune field was
flattened. By definition, the mean slope is given by the RDP (blue line, Eq.
5). Light and dark gray lines show the transport time scale for constant sand
fluxes of 10 and 30 ${\rm m}^{2}\,{\rm yr}^{-1}$, respectively. Highlighted
areas show the spring periods from March 21 to June 21.

###  6.4 Topographic surveys during experimental dune growth

From April 2014 to November 2017, we performed a series of topographic surveys
of the dunes developing from the flat sand bed using a ground-based laser
scanner Leica Scanstation C10 (Fig. S9a). To compare these different
measurements, a reference system of concrete posts was installed over the
entire experimental dune field (Figs. S9b-d). For each survey, the zone under
investigation is scanned from four different view points. Over the 42 months
of the experiment, the density of points varies from 472 to 2368 points/m2.

To map surface elevation, we select only elevation points within the 4-sided
polygon determined by reference points t6, t7, t8 and t9 (Figs. S9c-d). Within
this polygon, to avoid disturbances from the surrounding bedforms, we have
chosen a central rectangular area with a width of 48 m and a length of 82 m
(ABCD in Figs. S9c-d). The long side of this rectangle is oriented Northwest-
Southeast to align perpendicular to the final orientation of the dunes at the
end of the experiment. We remove the mean slope of this rectangular area by
adjusting a plane to the elevation data. This plane has always a southwest-
facing slope during the entire duration of the experiment. The residual
topography is shown for different times in Fig. 5A of the main manuscript.
Within the rectangular area, we study 2D transects parallel to the main axis
using a spacing of 1.4 m between two transects (see for example transect
$aa^{\prime}$ in Figs. S9c-d). For each of the 34 transects, we select the
elevation data points in a 0.2 m wide band on either side. All these points
are horizontally projected on their respective transect line. After
substracting the average slope and the mean elevation, we resample these data
to a regular spacing of 0.1, 0.25 and 0.35 m. Thus, we can use the fast
Fourier transform method to explore different values within the frequency
domain.

Fig. 5B of the main manuscript shows the elevation profiles with respect to
time for the transect $aa^{\prime}$ shown in Figs. S9c-d. Over the 42 months
of the experiment, the amplitude of the dunes varies from a few centimeters to
a few meters.

Figure S9:  Topographic surveys during experimental dune growth. (a) Ground
laser scanning in the landscape-scale experiment. (b) Ground view of three
posts of the local reference system. (c,d) Reference system (dots) and
selected area (blue rectangle) for the entire duration of the experiment. All
2D elevation profiles are taken parallel to the transect $aa^{\prime}$ (red
line). They are separated by a distance of 1.4 m. The elevation profiles of
the same transects at different times are used to study the dune instability.

###  6.5 Spectral analysis, amplitude and phase

For each transect and all topographic surveys, we use a fast Fourier transform
method to analyze the signal in the frequency domain (Figs. S10a,b). By
selecting only one frequency value, the individual contribution of each
wavelength is computed by an inverse Fourier transform (Fig. S10c). Thus, we
obtain a sine wave

$H_{\lambda_{i}}^{j}(t_{n})=h_{\lambda_{i}}^{j}(t_{n})\cos(k_{i}x+\phi^{j}_{\lambda_{i}}(t_{n}))$
(22)

for each wavelength $\lambda_{i}=2\pi/k_{i}$ of each transect $j$ at the
different times $t_{n}$ of the $n^{\rm th}$ topographic survey. Considering
the 34 different transects, the amplitude $h$ and the phase $\phi$ can be used
to estimate growth rates $\sigma_{\lambda_{i}}$ and phase velocities
$c_{\lambda_{i}}$ with respect to time, respectively. We perform this analysis
on the same elevation profiles sampled at different rates (0.1, 0.25 and 0.35
m) to verify that it has no effect on the spectral behavior.

Taking as an example transect No 20 in the middle of the selected area, Fig.
S11 shows the logarithm of the amplitude $A_{\lambda}^{20}$ of individual
wavelengths with respect to the dimensionless time scale. A sudden change in
behavior is observed between $t_{a}=10.7$ and $t_{a}=11.3$, i.e. from October
30 to November 12, 2014. Over this time period, there is an abrupt increase in
growth rate and then a more irregular behavior than during the initial phase.
As shown in Fig. 5C of the main manuscript, it is also during this period that
the mean amplitude of the bedforms begins to increase more rapidly. These
behaviors are interpreted as the transition between linear and non-linear
growth phases. This interpretation is supported by the observation of the
first slip faces downwind of the crest under the prevailing wind.

In order to quantify dune growth in the linear phase according to the dune
instability (see Sec. 6.1 and Eq. 20), we perform an exponential fit to the
amplitude data from $t_{a}=0$ to $t_{a}=10.7$, i.e. from April 10 to October
30, 2014, for each wavelength and each transect (see red lines in Fig. S11).
The agreement between the exponential regime and the data as well as the
variation of the exponential rate at different wavelengths indicate that there
is a coherent behavior at different length scales, which can be analyzed
thanks to a dispersion diagram.

Figure S10:  Spectral analysis of elevation profiles with respect to time. (a)
Two consecutive elevation profiles of transect No 18 during dune growth. (b)
Power spectral density of the elevation profiles using a fast Fourier
transform. Dashed lines show wavelengths $\lambda_{\\{1,2,3\\}}$ of
$\\{51.2,\,20.5,\,7.9\\}$ m. (c) The reconstructed elevation profiles using
individual wavelengths $\lambda_{\\{1,2,3\\}}$ and the inverse Fourier
transform. By definition, they take the form of sine waves
$H_{\lambda_{i}}=h_{\lambda_{i}}\cos(k_{i}x+\phi_{\lambda_{i}})$ where
$k_{i}=2\pi/\lambda_{i}$ is the wave number and $\phi_{\lambda_{i}}$ the
phase. Values of $h_{\lambda_{i}}$ and $\phi_{\lambda_{i}}$ at different times
can be used to estimate growth rate $\sigma_{\lambda_{i}}$ and phase velocity
$c_{\lambda_{i}}$, respectively.

###  6.6 Dispersion diagram of the growth rate

Fig. S11 shows the dispersion relation for transect No20 during the linear
phase from $t_{a}=0$ to $t_{a}=10.7$, i,e. from April 10 to October 30, 2014.
In this figure, the exponential rates of growth or decay of the different
wavelengths $\lambda_{i}$ are plotted with respect to the wave number
$k_{i}=2\pi/\lambda_{i}$. From the longest wavelengths (i.e. smallest wave
numbers), there is an increase in the growth rate. The maximum growth rate is
reached for $k_{\rm max}\approx 0.43\leavevmode\nobreak\ {\rm m}^{-1}$
($\lambda_{\rm max}\approx 14.6\leavevmode\nobreak\ {\rm m}$). For shorter
wavelengths, the growth rate is decreasing and the neutral mode is reached for
a value of $k_{0}$ of about $0.67\leavevmode\nobreak\ {\rm m}^{-1}$
($\lambda_{\rm c}\approx 9.3\leavevmode\nobreak\ {\rm m}$). Then, shorter
wavelengths have negative growth rates. These decay rates are associated with
small amplitudes of less than one centimeter, so that their rapid variations
can not be captured given the resolution of our topographic data.

Fig. S12a shows the dispersion relations of the growth rate for all transects
and the three sampling rates. The mean and standard deviation of these growth
rates are plotted with respect to the wave number in Fig. S12b. Using these
data, we find $k_{\rm max}\approx 0.45\leavevmode\nobreak\ {\rm m}^{-1}$
($\lambda_{\rm max}\approx 14\leavevmode\nobreak\ {\rm m}$) and $k_{\rm
c}\approx 0.7\leavevmode\nobreak\ {\rm m}^{-1}$ ($\lambda_{\rm c}\approx
9\leavevmode\nobreak\ {\rm m}$).

Theoretically, the growth rate writes

$\sigma(k)=Q\,k^{2}\,\dfrac{B-A\,k\,l_{\rm sat}}{1+\left(k\,l_{\rm
sat}\right)^{2}}$ (23)

where $A$ and $B$ are the aerodynamic parameters (Sec. 5), $l_{\rm sat}$ the
saturation length (Sec. 4), $Q$ a constant flux proportional to $u_{*}^{2}$.
This equation does not take into account the transport threshold and is
therefore theoretically valid only within the limit of strong winds when the
sand flux is large. However, the role of the transport threshold is included
in our estimates thanks to the dimensionless time scale and a prefactor
correcting for the intensity of the sand flux. In practice, the growth rate
derived from the elevation profiles scales as $\sigma(k)/Q$ and we consider a
prefactor $\langle u_{*}\rangle^{2}/(\langle u_{*}\rangle^{2}-u_{\rm th}^{2})$
that integrates all the variations of wind strength. Thus, the parameter $A$
and $B$ in Eq. 18 and in Eq. 23 are compatible with one another.

The dispersion relation using Eq. 23 and the values of $\\{l_{\rm
sat},\,A,\,B\\}=\\{0.95,\,3,\,1.5\\}$ measured independently in the field are
plotted in orange in Figs S11 and S12 as in Fig. 6B of the main manuscript.
The blue curve in these figures corresponds to a dispersion relation with the
best-fit values of $\\{l_{\rm sat},\,A,\,B\\}=\\{0.7,\,1.96,\,0.96\\}$. In
addition, to investigate the sensitivity of the dispersion relation to the
values of $\\{l_{\rm sat},\,A,\,B\\}$, the shaded area in Fig. 3c of the main
manuscript shows the best-fit values of $A$ and $B$ using values of $l_{\rm
sat}\in[0.5;\,1.1]$ in Eq. 23. All these results underline the consistency of
the experimental dispersion diagram derived from the topographic data not only
with the theory but also with our independent measurements of flow and
transport properties in the field.

Figure S11:  Growth rate of individual wavelengths. Considering the elevation
profiles of transect No20 and a sampling rate of 25 cm, figures show the
logarithm of the amplitude $h_{\lambda}$ of individual wavelengths with
respect to the dimensionless time scale; $t_{a}=0$ is set for April 10, 2014.
Two different regimes are observed before and after $t_{a}=10.7$, November 30,
2014. There are associated with the linear and non-linear growth phases. The
growth rate of each wavelength is determined by an exponential fit performed
during the linear phase ($0\leq t_{a}\leq 10.7$, dashed lines). As a
synthesis, the bottom right figure shows the dispersion relation, i.e. the
growth rate with respect to the wave number $k=2\pi/\lambda$. Solid lines are
the same dispersion relations as in Fig. 6B of the main manuscript: $\\{l_{\rm
sat},\,A,\,B\\}=\\{0.95,\,3,\,1.5\\}$ (orange) and $\\{0.7,\,1.96,\,0.96\\}$
(blue). Note that small wavelengths have amplitude of less than one
centimeter. These amplitudes are comparable to the noise level of our
topographic measurements, which prevents the estimation of the exponential
decrease in growth rate for short wavelengths ($\leq$ 8 m).

Figure S12:  Dispersion diagram of the growth rate for all transects at
different sampling rates. (a) Growth rates during the linear phase, i.e. from
April 10, 2014 to November 30, 2014, with respect to the wave number
$k=2\pi/\lambda$ for each transect and sampling rates of 0.1 m (blue), 0.25 m
(yellow) and 0.35 m (green). (b) Mean and standard deviation of the growth
rates averaged on all transects. Solid lines are the same dispersion relations
as in Fig. 6B of the main manuscript and in Fig. S11: $\\{l_{\rm
sat},\,A,\,B\\}=\\{0.95,\,3,\,1.5\\}$ (orange) and $\\{0.7,\,1.96,\,0.96\\}$
(blue).

###  6.7 Dispersion diagram of the phase velocity

During the linear-growth phase, the analysis of the variation of the phase
$\phi(\lambda_{i})$ (Eq. 22) as a function of time do not at present provide
conclusive evidence about the dispersion diagram of the phase velocity. An
obvious limitation comes from the time delay between two consecutive
topographic surveys. Indeed, the phase shift can be too large to be accurately
computed, even by folding the data modulo the period from large to short
wavelengths. The orientation and magnitude of the resultant sand flux
associated with each time interval is another obvious issue. In fact, wind
reversals and the subsequent back and forth of the incipient bedforms were
frequent during the experiment. Variation of the phase data during periods of
strong winds and small resultant sand flux are particularly difficult to be
interpreted.

To simplify the problem, we can also estimate dune migration distance by cross
correlation between two elevation profiles over a time interval over which the
wind is almost unidirectional, from April 25 to May 13, 2014, (Figs. S13a,b).
The dune migration rate can then be computed using the dimensionless time
scale to be compared to theoretical dispersion relations. With the same
notation as in Eq. 23, dispersion relations write

$c(k)=Q\,k\,\dfrac{A+B\,k\,l_{\rm sat}}{1+\left(k\,l_{\rm sat}\right)^{2}}.$
(24)

Considering the dimensionless time scale and the same prefactor as for the
growth rate, Fig. S13c shows the dispersion relations of the phase velocity
using the same parameters as in Figs. S11 and S12. It also shows the phase
velocity of the most unstable wavelength $\lambda_{\rm max}$, the mean value
of the migration rate averaged over all transects and the corresponding
standard deviation. As expected, the migration rate derived from the cross
correlation is close to the phase velocity of the most unstable wavelength.

Figure S13:  Dispersion diagram of the phase velocity. (a) Two consecutive
elevation profiles of transect No 16 during dune growth. (b) The maximum value
of the cross correlation diagram between the two elevation profiles determines
the dune migration distance $\Delta l=3.25\leavevmode\nobreak\ {\rm m}$ from
April 25 to May 13, 2014 (i.e. $\Delta t_{a}=1.77$). (c) Dispersion diagram of
the phase velocity using the same parameters as in Fig. 6B of the main
manuscript and in Figs. S11 and S12b: $\\{l_{\rm
sat},\,A,\,B\\}=\\{0.95,\,3,\,1.5\\}$ (orange) and $\\{0.7,\,1.96,\,0.96\\}$
(blue). Vertical dashed lines show the most unstable wavelength $\lambda_{\rm
max}$ in both cases. Using cross correlation, the horizontal dashed lines
shows the mean migration velocity $\Delta l/\Delta t_{a}$ averaged over all
transects. The shaded area has a width of two standard deviation using the
same data. There is a quantitative agreement between $c(\lambda_{\rm max})$
and the migration rate derived from $\Delta l/\Delta t_{a}$. This is possible
because during this period of time of the linear growth phase the wind regime
is almost unidirectional (see inset).

## References

  * [1] Iversen, J. D. & Rasmussen, K. R.  The effect of wind speed and bed slope on sand transport.  _Sedimentology_ 46, 723–731 (1999).
  * [2] Ungar, J. & Haff, P.  Steady state saltation in air.  _Sedimentology_ 34, 289–299 (1987).
  * [3] Durán, O., Claudin, P. & Andreotti, B.  On aeolian transport: Grain-scale interactions, dynamical mechanisms and scaling laws.  _Aeolian Research_ 3, 243–270 (2011).
  * [4] Fryberger, S. G. & Dean, G.  Dune forms and wind regime.  _A study of global sand seas_ 1052, 137–169 (1979).
  * [5] Pearce, K. I. & Walker, I. J.  Frequency and magnitude biases in the ‘Fryberger’ model, with implications for characterizing geomorphically effective winds.  _Geomorphology_ 68, 39–55 (2005).
  * [6] Tsoar, H.  Sand dunes mobility and stability in relation to climate.  _Physica A_ 357, 50–56 (2005).
  * [7] Bagnold, R.  The physics of wind blown sand and desert dunes: London (1941).
  * [8] Sauermann, G., Kroy, K. & Herrmann, H. J.  Continuum saltation model for sand dunes.  _Physical Review E_ 6403, art. no.–031305 (2001).
  * [9] Narteau, C., Zhang, D., Rozier, O. & Claudin, P.  Setting the length and time scales of a cellular automaton dune model from the analysis of superimposed bed forms.  _J. Geophys. Res._ 114 (2009).
  * [10] Andreotti, B., Claudin, P. & Douady, S.  Selection of dune shapes and velocities: Part 1: Dynamics of sand, wind and barchans.  _Eur. Phys. J. B_ 28, 321–339 (2002).
  * [11] Charru, F.  Selection of the ripple length on a granular bed.  _Phys. Fluids_ 18, 121–508 (2006).
  * [12] Elbelrhiti, H., Claudin, P. & Andreotti, B.  Field evidence for surface-wave-induced instability of sand dunes.  _Nature_ 437, 720–723 (2005).
  * [13] Fourrière, A., Claudin, P. & Andreotti, B.  Bedforms in a turbulent stream: formation of ripples by primary linear instability and of dunes by nonlinear pattern coarsening.  _Journal of Fluid Mechanics_ 649, 287–328 (2010).
  * [14] Hersen, P., Douady, S. & Andreotti, B.  Relevant length scale of barchan dunes.  _Physical Review Letters_ 89 (2002).
  * [15] Andreotti, B.  The song of dunes as a wave-particle mode locking.  _Physical Review Letters_ 93 (2004).
  * [16] Andreotti, B., Claudin, P. & Pouliquen, O.  Measurements of the aeolian sand transport saturation length.  _Geomorphology_ 123, 343–348 (2010).
  * [17] Pähtz, T., Kok, J. F., Parteli, E. J. & Herrmann, H. J.  Flux saturation length of sediment transport.  _Physical review letters_ 111, 218002 (2013).
  * [18] Jackson, P. & Hunt, J.  Turbulent wind flow over a low hill.  _Q. J. R. Meteorol. Soc._ 101, 929–955 (1975).
  * [19] Sykes, R. I.  An asymptotic theory of incompressible turbulent boundary layer flow over a small hump.  _Journal of Fluid Mechanics_ 101, 647–670 (1980).
  * [20] Claudin, P., Wiggs, G. & Andreotti, B.  Field evidence for the upwind velocity shift at the crest of low dunes.  _Boundary-layer meteorology_ 148, 195–206 (2013).