# Emergence of wind ripples controlled by mechanics of grain-bed impacts

Conner W. Lester, A B Murray, Orencio Duran, Bruno Andreotti, Philippe Claudin

► To cite this version:

Conner W. Lester, A B Murray, Orencio Duran, Bruno Andreotti, Philippe Claudin. Emergence of wind ripples controlled by mechanics of grain-bed impacts. Nature Geoscience, 2025, 18 (4), pp.344-350. 10.1038/s41561-025-01672-w. hal-05053701

HAL Id: hal-05053701

https://hal.science/hal-05053701v1

Submitted on 2 May 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Emergence of wind ripples controlled by mechanics of grain-bed impacts

C. W. Lester,¹,* A. B. Murray,¹ Orencio Duran,² B. Andreotti,³ and P. Claudin⁴

¹Earth and Climate Sciences, Duke University, Durham, NC 27708, USA

²Department of Ocean Engineering, Texas A&M University, College Station, TX 77843-3136, USA

³Laboratoire de Physique de l'Ecole Normale Supérieure, PSL, Sorbonne Université,
Université Paris Cité, CNRS, 24 rue Lhomond 75005 Paris, France

⁴Physique et Mécanique des Milieux Hétérogènes, PMMH UMR 7636 CNRS,
ESPCI Paris, PSL Research University, Sorbonne Université,
Université Paris Cité, 7 quai St Bernard, 75005 Paris, France

ABSTRACT

Periodic sediment patterns have been observed on Earth in riverbeds and sand and snow deserts, but also in other planetary environments. One of the most ubiquitous patterns, familiar wind or ‘impact’ ripples, adorns sand beaches and arid regions on Earth. The observation of aeolian impact ripples on Mars the same size as their terrestrial counterparts despite a thinner atmosphere raises questions about their formation. Here we show in a numerical simulation that the emergent wavelength of impact ripples is controlled by the mechanics of grain-bed impacts and not the characteristic trajectories of grains above the bed. We find that the distribution of grain trajectories in transport is essentially scale-free, invoking the proximity of a critical point and precluding a transport-related length scale that selects ripple wavelengths. By contrast, when a grain strikes the bed, the process leading to grain ejections introduces a collective granular length scale that determines the scale of the ripples. We propose a theoretical model that predicts a relatively constant ripple size for most planetary conditions. In addition, our model predicts that for high-density atmospheres, such as on Venus, or for sufficiently large sand grains on Earth, impact ripples propagate upwind. Although wind-tunnel and field experiments are needed to confirm the existence of such ‘antiripples’, we suggest that our quantitative model of wind-blown sediment transport may be used to deduce geological and environmental conditions on other planets from the sizes and propagation speeds of impact ripples.

THIS IS A PREPRINT; IT HAS NOT BEEN PEER REVIEWED BY A JOURNAL.

A version of this preprint was published at Nature Geoscience on April 8th, 2025. See https://doi.org/10.1038/s41561-025-01672-w.

<footnote>* Corresponding author contact: conner.lester@duke.edu</footnote>
Many approximately periodic sedimentary patterns result from the interaction between bedform and flow, via modulation of transport [1, 2]. There is already a long scientific history of relating bedforms observed in different environments to physical parameters such as fluid density, viscosity and velocity, and grain density and size [1–4]. The search for analogs between different planets involves correctly identifying the dynamical mechanisms. For example, the meter scale dunes on Mars (Fig. 1b) are not the counterparts of centimeter scale wind ripples on Earth (Fig. 1a), but are more similar to ripples in water and sand dunes more broadly [4]. While dunes, water ripples and meter-scale Martian bedforms arise because fluid forces on the bed vary across the bedform, wind ripples, or ‘impact ripples’, form under conditions in which grains at the bed surface are set in motion by high energy impacts of other grains in transport—as opposed to directly impelled by the fluid. On both Earth and Mars, impact ripples are centimeter scale (Fig. 1a,b). However, experimental results suggest that under some conditions wind modulation in response to bed undulations can lead to the coevolution of impact and hydrodynamic ripples of similar scale [5]. To disentangle these effects we choose to analyze impact-ripple formation in the absence of wind modulations.

The current understanding of impact ripples often attributes their size (wavelength) to the characteristic hoplength of sand grains in transport [3, 6–9], thought to be the only relevant length scale. However, this explanation implies that impact ripple wavelengths should tend to be larger on planets where lower density atmospheres lead to larger hoplengths, such as Mars, and smaller for thicker atmospheres with smaller hoplengths, such as Venus—in contrast with observations [4, 10–13]. We investigate ripple emergence using a discrete element model (DEM) with a two way coupling between grains and hydrodynamics [14], able to reproduce observed transport characteristics and impact ripples over a wide range of atmospheric environments [5, 9, 15, 16] (see Supplemental Videos). Simulations show that the hoplength distribution is essentially scale-free over nearly the whole range of hoplengths $\ell$ (Fig. 1d; [9]). This implies that no characteristic hoplength scale relevant for ripple wavelength selection exists [6], and that our understanding of ripple formation and aeolian sediment transport is flawed. In this paper we explore the consequences of the scale-free nature of sediment transport and propose a new dynamic length scale resulting from grain-bed impacts that determines impact ripple wavelengths across planetary environments.

## Criticality Driven Surface Instability

Figure 1d shows that grain trajectories vary continuously between the limits of very long, high-energy trajectories and more numerous very short hops, in contrast with the standard conceptual model of wind-driven transport that divides grains in motion into two distinct populations, ‘saltons’ and ‘reptons’ [6–8]. Saltons represent long trajectory, high energy grains that impact the sand surface and eject low energy reptons, which make relatively short hops. Counter to the bimodal salton-repton picture, the steady-state distribution of hoplengths, $P_0(\ell)$, measured from the DEM can be approximated as a $\sim \ell^{-1}$ power-law with an exponential cutoff for exceedingly long hoplengths set by grains that make it outside of the transport layer and whose maximum hoplengths are set by the wind speed (Fig. 1e). Wind tunnel observations [17, 18] also report a roughly $P_0 \sim \ell^{-1}$ scaling with an exponential tail (Fig. 1e inset).
<figure><img src="image_2.png" /><figcaption>FIG. 2. Trajectory-ripple interactions lead to flux modulations, and impact-ejection lag phase shifts the modulated flux. (a) DEM trajectories over a rippled surface showing only hoplengths $l$ around one wavelength $\lambda$—illustrating the modulation of the impact flux via trajectory amplification. Panels (b) show the processes that lead to trajectory and resulting flux modulations. (b1) The mechanism for amplification and focusing of trajectories around the upwind inflection point due to geometrical modulations of hoplengths $\Delta l(l, x, t)$. Red amendments to trajectories show positive or negative changes to hoplengths $l$ based on where the grain impacts for a given take off location. Note that the solid red trajectory of $l \sim \lambda$ is the attractive in the sense that, with iterations, grains will tend towards hoplengths of $l \sim \lambda$. Contrary, the dashed trajectory is repulsive with iterations. Panel (b2) shows how the impact flux modulation is achieved by successive iterations of trajectories starting from a uniform impact flux $\phi_0$ (note the increased flux at the upwind inflection point and decreased flux at the downwind inflection point). The final panel (b3) shows how the impact flux is modulated on a per hoplength basis (the impact flux density). The red line $l = \lambda$ reflects the impact flux density of trajectories in panel (a). The sum over all impacting hoplengths results in the total impact flux modulation as in the solid line in (b2). Panels (b) are produced for $\lambda = 200d$ as in (a) which give rise to the specific phase shifts observed (see Supplementary Information). Panels (c) shows how the granular dynamics of impact lead to a granular length scale that can lead to a phase shift in the vertical flux. (c1) An impact from the DEM of a high energy grain (red) with the surface. $\epsilon = x_\uparrow - x_\downarrow$ is the location of an ejection relative to the impact location. Panel (c2) shows DEM data for the number density of ejections and rebounds per impact of a given hoplength $n_\epsilon(\epsilon|l)$ versus rescaled $\epsilon/d$ for various ranges of impacting hoplengths (a proxy for impact energy). The average impact lag distance $\langle\epsilon\rangle$ has the main effect of simply phase shifting the ejection flux downwind—as illustrated in the schematic of (c3).</figcaption></figure>

collisions occur on a bed of grains on the verge of motion. Thus the power-law distribution of hoplengths may result from the collisional process on an extremely sensitive grain arrangement [22]. (Power-law distributions of temporal fluctuations in the number of grains in transport also suggest critical behavior [20]; Methods Fig. M1.)

Because grains impact the surface at low angles, sediment transport over ripples tends to concentrate on upwind faces, leading to an increase in ejections from the face that get distributed across the ripples crest (Fig. 1d, 2a,b). This constitutes the destabilizing mechanism that makes the ripple grow [6–8]. Additionally, the density of impacts on upwind faces tends to be iteratively amplified moving downwind (Fig. 2b). The amplification occurs for a subpopulation of grains leaving an upwind face with hoplengths close to one ripple wavelength (Fig. 2a): Hoplengths just longer than a wavelength are shortened, and shorter hoplengths are extended, depending on the relative difference in elevation between lift-off and landing (Fig. 2b; Methods). This focusing leads to a further concentration of impacts
<figure><img src="image_3.png" /><figcaption>FIG. 3. Scale-free properties of sediment transport. The DEM suggests that the majority of trajectories ($\ell \lesssim 10\langle\ell\rangle$) can be characterized by power-law functions for the hoplength distribution $P_0(\ell)$, the replacement capacity $n_\uparrow(\ell)$ (the number of ejections and rebounds per impact) and the mean impact-ejection lag distance $\bar{\epsilon}(\ell)$—three among many other dynamical functions that show power-law behavior in the DEM. The colors show DEM data corresponding to conditions adjusted to represent Mars, Earth, Titan and Venus by varying the grain-fluid density ratio $\rho_p/\rho_f = (10^5, 2000, 200, 50)$, respectively. All simulations are done in rarefied transport conditions with $d \simeq 100\mu m$ (Methods). Inset triangles in each plot are power law references. We note that the average lag distance $\int \bar{\epsilon}(\ell)P_0(\ell)d\ell \sim 2d$ for all four conditions.</figcaption></figure>

on the subsequent upwind face, and therefore to a further increase in the number of grains rebounding and ejected from there, completing the feedback.

As a consequence of the scale-free hoplength distribution, the amplification of impacts on upwind slopes occurs for any wavelength perturbation, and the weighting toward small hops tends to make the smallest wavelengths grow fastest (see Supplementary Information, or SI, Section II). The fact that ripples in nature emerge at finite characteristic wavelengths (much greater than the grain size $d$) suggests that some mechanism stabilizes the smallest wavelengths.

## Granular Dynamics Wavelength Selection

Analysis of DEM results reveal a spatial lag, $\epsilon = x_\uparrow - x_\downarrow$, between the location of an impact $x_\downarrow$ and a resulting ejection at $x_\uparrow$ (Fig. 2c). The momentum from an impact is distributed through a complex network of dynamic force chains [23] around the impact location, and because grains impact the bed at low angles this momentum is preferentially directed downwind. Thus, grains are often ejected several grain diameters away from an impact site— resulting in a length scale $\epsilon$ that arises from collective granular dynamics. The spatial lag $\epsilon$ provides the scale that determines the wavelength below which ripples cannot grow, ultimately selecting the emergent ripple size.

The lag in the locations of ejections relative to impacts shifts the position of the erosion-deposition transition point relative to the bedform—as the “saturation length” does for hydrodynamic bedforms [2]. When the perturbation wavelength is sufficiently long, the lag does not play a significant role. However, when the wavelength is small enough, the lag shifts the transition point downwind of the crest, so that the crest erodes and the bedform decays. The fastest growing wavelength emerges from the competition between the destabilizing mechanism, which favors the growth of small wavelengths, and the stabilizing mechanism arising from the lag, with the result that the emergent wavelength scales with the average impact-ejection lag distance $\langle\epsilon\rangle$ (Methods).
Analytical Model

In the Methods and Supplementary Information (SI) we derive a generalized model for ripple growth which couples modulated grain trajectories—and resulting modulated sediment fluxes (Fig. 2b)—with the nonlocal impact-ejection process (Fig. 2c). To determine what selects initial ripple wavelengths and propagation speeds we perform a linear stability analysis to derive an expression for the ripple growth rate $\sigma(k)$ and propagation speed $c(k)$ as functions of perturbation mode $k = 2\pi/\lambda$, where $\lambda$ is the perturbation wavelength (the ‘dispersion relation’ for ripple growth). This analysis reveals that ripple growth is a function of some key transport quantities measured on a flatbed: the hoplength distribution $P_0(\ell)$, the average number of ejections and rebounds per impact, or the replacement capacity $n_\uparrow(\ell)$, and the average lag distance per impact $\bar{\varepsilon}(\ell)$, which all show scale-free/power-law behavior in the DEM for the majority of trajectories (Fig. 3), again hinting at the criticality of sediment transport [9, 22].

In order to measure the replacement capacity $n_\uparrow(\ell)$ in the DEM *during transport*, we link individual grains leaving the sand bed (crossing a horizontal plane at $z = 0.7d$ from the surface) to grains that have impacted that area ($dx = 20d$) in a time window before the crossing ($dt = 2\sqrt{d/g}$). Values for the crossing elevation and the space and time windows are chosen to satisfy the constraint that most ejecta are paired with an impact. If the particle leaving the surface has impacted the surface within that time and space window then we count it as a rebound. Otherwise, we search for the highest energy impactor that likely gave rise to that ejecta. And $\varepsilon = x_\uparrow - x_\downarrow$ is simply the measured displacement between impact and ejection/rebound.

Using the DEM measurements of flatbed quantities in Figure 3, we find the derived dispersion relation agrees with numerical impact ripples that emerge in the DEM (Fig. 4, markers), and is broadly consistent with ripple observations on Earth, as well as other planets (Fig. 4, shaded regions). The numerical ripples produced in the DEM have been directly compared to wind-tunnel experiments for a wide range of atmospheres (from Earth to Mars-like as well as Venus-like conditions) [5, 13]. Therefore, the fact that the analytical model can reproduce the numerical ripple emergence suggests that the mechanisms accounted for in this new conceptual and analytical model capture the main physics of the system. This agreement also implies that ripple emergence is tied to the critical nature of sediment transport, because the dispersion relations of Figure 4 (black lines) are derived using the power-law scalings observed during transport (Fig. 3; Methods).

In this new analytical framework, ripple wavelengths and speeds arise from the modulation of the impacting flux of grains, which strongly depends on the hoplength $\ell$ due to the focusing effect (Fig. 2b). This is encoded in two dimensionless quantities $\mathcal{R}(\ell)$ and $\mathcal{S}(\ell)$ which are rescaled components of the modulated impact flux density of hoplengths in-phase and in-quadrature with respect to the topography (Methods). The in-phase modulation $\mathcal{R}(\ell)$ quantifies the contribution of trajectories of hoplength $\ell$ to the amplification of the impacting flux on ripple crests, and the out-of-phase modulation $\mathcal{S}(\ell)$ quantifies the contribution on ripple faces (e.g., Fig. 2b). Because these modulations result from topography-induced focusing effects (Fig. 2b), $\mathcal{R}(\ell) \simeq \mathcal{R}(\ell/\lambda)$ and $\mathcal{S}(\ell) \simeq \mathcal{S}(\ell/\lambda)$ primarily depend on the ratio $\ell/\lambda$ of their hoplength to the wavelength. As a first approximation $\mathcal{R} \sim \sin k\ell$ and $\mathcal{S} \sim 1 + \cos k\ell$ (Methods
<figure><img src="image_4.png" /><figcaption>FIG. 4. Ripple growth rate and speed curves for Mars (red), Earth (blue), Titan (green) and Venus (gold). Circle markers are DEM experiments starting with initial bed modulations of a given wavelength $\lambda$ (see Methods). Filled circles are positive values and open circles represent negative values. Diamond markers are DEM measurements from ripple evolution from a flat initial surface condition (see space-time plots in (b)). Black lines are predictions from the new analytical model (Methods and SI) using power law approximations for input functions (Fig.'s 3, M2). Dashed black lines are negative values. (a) Ripple growth rates $\sigma$ rescaled by the homogeneous vertical flux $\phi_0$ and the grain size $d$ versus rescaled wavelength $\lambda/d$ predict fastest growing wavelengths of roughly $\lambda \sim 200d$ for all conditions. (b) Rescaled ripple speeds $c/\phi_0$ versus wavelength $\lambda$ rescaled by the average hoplength $\langle l \rangle$. Inset plots in (b) show space-time evolution of DEM ripple simulations. The grey lines for Earth in (a,b) show predictions from the salton-repton model [6], using the power-law like distribution of hoplengths $P_0(l)$ (SI Section IID). Shaded backgrounds in all plots show ranges of observed ripple sizes $\lambda$ from wind tunnel and field experiments (Mars [4, 5], Earth [5, 11, 12, 24–26], Venus [13]) which are in broad agreement with predictions. We note that Venus wind tunnel experiments [13] did not report ripple propagation speeds—even though they were likely present—thus we can not confirm if they observed antiripples $c < 0$, as predicted by our model. Additionally there is no observational data for Titan (wind tunnel or field/satellite) thus the vertical dashed line, $\lambda \sim 200d$ and $c \sim 0$, represents predictions from our model.</figcaption></figure>

Fig. M3), which are maximal for $l \simeq \lambda/4$ (representing grains landing on the crest after being ejected on the upwind flank) and $l \simeq \lambda$ (representing amplified trajectories traveling from one upwind flank to the next) respectively. Our analytical model (Methods) predicts that fastest growing wavelength and associated propagation speed scale as

$$ \lambda \sim \frac{\langle n_{\uparrow} S \rangle}{\langle (1 - n_{\uparrow}) R \rangle} \langle \varepsilon \rangle \quad \text{and} \quad c \sim \langle (n_{\uparrow} - 1) S \rangle \phi_0, \qquad (1) $$

which shows that the wavelength scale is set by the ratio of the stabilizing mechanism, which scales with the average number of grains leaving from the ripple face $\langle n_{\uparrow} S \rangle > 0$ and the shift $\langle \varepsilon \rangle$, and the destabilizing mechanism, which scales with the average number of grains deposited on the crest $\langle (1 - n_{\uparrow}) R \rangle > 0$ (Methods). The propagation speed scales with the homogeneous vertical flux at the surface $\phi_0$ and the number of grains eroded from (or deposited on) the ripple faces $\langle (n_{\uparrow} - 1) S \rangle$, which may be positive or negative depending on the ratio $\langle l \rangle/\lambda$ (Fig. 4; Methods Fig. M4),
where $\langle \ell \rangle$ is the average hoplength set by the exponential cutoff of $P_0(\ell)$ (Fig. 3).

## Antiripples

Negative propagation speeds represent the novel prediction that under certain conditions impact ripples will propagate upwind as antiripples. Upwind propagation will occur if there is net deposition on the upwind inflection point (Eq. 1). When the hoplength scale is of the order of the ripple wavelength ($\langle \ell \rangle \sim \lambda$) there are many modulated high-energy trajectories (i.e., $\ell \sim \lambda$) that cause erosion of the upwind face, $\langle n_\uparrow S \rangle > \langle S \rangle$, and downwind propagation ($c > 0$, Eq. 1). In contrast, in the limit $\langle \ell \rangle \ll \lambda$, most grains have hoplengths $\ell$ smaller than $\lambda$. Therefore on upwind slopes nearly all trajectories originate locally, and have their hoplengths shortened (Fig. 2b), tending to cause a reduction in the number of ejections per impact, $n_\uparrow(\ell)$ (Fig. 3). Thus, more grains land around the upwind inflection point than leave, $\langle S \rangle > \langle n_\uparrow S \rangle$, which results in net accretion on the upwind face and upwind propagation as antiripples ($c < 0$, Eq. 1; Methods Fig. M4).

Therefore, our analysis predicts that antiripples will emerge under conditions in which the fastest growing wavelength, which is ultimately selected by the granular length scale $\langle \epsilon \rangle$, is much larger than the hoplength scale $\langle \ell \rangle$. This is predicted to occur when the density ratio between the grains and the fluid $\rho_p/\rho_f$ is small and transport is tending to transition to bedload (smaller hoplengths), as on some other planetary bodies such as Titan and Venus (Fig. 4b), and may also occur on Earth for gravel sized grains where the hoplengths are also sufficiently short.

## Implications

Impact ripples occur across the globe and solar system. In our new model for the emergence of these patterns, we find that rather than ripple wavelengths being selected by a sediment transport length scale, they are selected by a bed deformation length scale, $\langle \epsilon \rangle$. This length scale is analogous to the impact-ejection lag length in nanometer-scale ion sputtering (erosion) processes which selects the wavelengths of micron-scale ripples on surfaces bombarded by ions (at low temperatures) [27, 28]. Additionally, the mechanism that allows $\langle \epsilon \rangle$ to determine wavelengths that are orders of magnitude larger (via scale-free impact forcing) is analogous to the amplification of the sediment-flux saturation length $L_{sat}$ which selects dune sizes (scale-free fluid forcing [1, 2]).

The importance of the length scale $\epsilon$ provides insights into long standing problems of ripple sizes on other planetary bodies. Because this bed deformation length is a granular phenomenon dependent mostly on grain properties, such as $d$, and depending little on the structure of the fluid or atmospheric density, we may expect that ripple wavelengths on other planets with the same grain composition should have roughly the same impact ripple sizes—which is what we find when comparing Earth to Venus and Martian impact ripples (Fig. 4; [4, 5, 13]).

While wavelengths are selected by grain-bed interactions, ripple propagation speeds are almost uniquely determined by the vertical flux scale $\phi_0$—set by the wind speed—and the scale of the hoplength distribution $\langle \ell \rangle$—set by the
atmospheric density and grain density and size. This suggests that while wavelengths may not vary much from planet to planet, propagation speeds may vary significantly—even from downwind to upwind propagating *antiripples* (Fig. 4b). Predictions of negative speed antiripples on Earth, with sufficiently coarse or light grains, can be tested in wind tunnels and field experiments.

Subsequent work should address how wind modulation affects the formation of impact ripples, as well as the convolution of impact and hydrodynamically forced bedforms [5]. Additionally, pattern coarsening from nonlinear 3-dimensional dynamics over long time scales modify the ripples that initially emerge from the 2-dimensional linear dynamics we focus on here [29–31]. Although the density of defects in the initial pattern determines how much wavelength can increase [31, 32], wavelengths typically only coarsen by approximately a factor of two or so based on wind tunnel observations [10, 12] and the distribution of ripple wavelengths observed on Earth and Mars [4]. Thus it is unlikely that coarsening effects are significant enough to erase the memory of the linear instability. Additionally, previous modeling [9] and experimental [12] work indicates that the initially fastest growing wavelength can increase moderately with wind speed. Wavelengths increase when wind speed, sediment flux, and impact densities are high enough that a continuously excited ‘gaseous layer’ forms at the bed surface [9]. Impact ripples can form for wind speeds only slightly above the threshold for sediment transport, when impacts are rarefied enough that no gaseous layer exists. We have focused here on conditions near the threshold, to explore the most general, fundamental mechanisms causing ripple formation. Possible future work could derive parametrizations analogous to those in Fig. 2 for higher wind speed, enabling an investigation of the relationship between wavelength and wind speed using our analytical model (Methods and SI).

Future efforts to compare model predictions to observations will test the new model of ripple formation, but more importantly they will also implicitly test the implications that scale-free transport dynamics (Fig. 3) and granular length scales are important fundamental aspects of aeolian sediment transport, generally.
# Methods

## Surface Evolution

To quantitatively explore the conceptual model described in the main text, we consider the evolution of the surface elevation Z(x, t) using the conservation of mass, which can be alternately expressed in terms of gradients in horizontal flux or in terms of vertical fluxes (SI I):

$$ \partial_t Z(x, t) = \varphi_{\downarrow}(x, t) - \varphi_{\uparrow}(x, t), \qquad (M1) $$

where $\varphi_{\downarrow}(x, t)$ and $\varphi_{\uparrow}(x, t)$ are the vertical fluxes of grains respectively reaching the bed and leaving the bed, normalized by the bed packing fraction $\phi_b$.

## Linking Impacts to Ejections

The impact flux density $\psi(\ell, x_{\downarrow})$ is defined such that $\psi(\ell, x_{\downarrow})d\ell dx_{\downarrow}$ is the volume of the grains (packed at the bed volume fraction $\phi_b$) arriving per unit time in the interval $[x_{\downarrow}, x_{\downarrow} + dx_{\downarrow}]$ after a hop of length between $\ell$ and $\ell + d\ell$. It can be expressed as a function of the distribution of the hoplengths $\ell$ on a flatbed, noted $P_0(\ell)$ (Fig. 1e, 3), and of the hoplength modulation, noted $\Delta\ell(\ell, x_{\downarrow})$, which result from surface perturbations (Fig. 2b). From the conservation of mass, $\psi(\ell, x_{\downarrow})$ can be linked to upwind ejections/rebounds, $\varphi_{\uparrow}(x_{\downarrow} - \ell)$, by the expansion (SI I C)

$$ \psi(\ell, x_{\downarrow}) \sim \varphi_{\uparrow}(x_{\downarrow} - \ell)P_0(\ell) - \varphi_0(\partial_{\ell} + \partial_{x_{\downarrow}})\Delta\ell(\ell, x_{\downarrow})P_0(\ell) \qquad (M2) $$

where $\varphi_0$ is the vertical flux in the homogeneous and steady-state, for which $\varphi_{\downarrow}(x) = \varphi_{\uparrow}(x) = \varphi_0$ and $\psi(\ell, x) = \psi_0(\ell) = \varphi_0 P_0(\ell)$. Stated simply, the impact flux density of grains is the translation of the ejection flux of grains $\varphi_{\uparrow}(x_{\downarrow} - \ell)$ traveling a distance $\ell = x_{\downarrow} - x_{\uparrow}$ determined by the homogeneous distribution $P_0(\ell)$, plus the flux modulation associated with the changes in trajectories $\Delta\ell(\ell, x_{\downarrow})$. Integrating this equation over $\ell$ links the total impact flux $\varphi_{\downarrow}$ to the vertical flux $\varphi_{\uparrow}$ by adding up the flux contributions from all trajectories:

$$ \begin{aligned} \varphi_{\downarrow}(x_{\downarrow}) &= \int \psi(\ell, x_{\downarrow})d\ell && (M3) \\ &= \int \varphi_{\uparrow}(x_{\downarrow} - \ell)P_0(\ell)d\ell - \varphi_0 \partial_{x_{\downarrow}} \langle \Delta\ell(x_{\downarrow}) \rangle \end{aligned} $$

The second term involves the spatial gradient of the average change in hoplengths $\langle \Delta\ell(x_{\downarrow}) \rangle$ (SI II A).
## Hoplength Modulations

In the presence of a modulated surface, the hoplength $\ell$ is not in general the same as the flat bed hoplength $\ell_0$ and can be written as

$$ \ell = \ell_0 + \Delta\ell, \qquad (M4) $$

where hoplength changes $\Delta\ell(\ell, x, t)$ are a function of trajectory interactions with surface geometry (at lift-off $x_\uparrow$ and on impact $x_\downarrow$) and wind modulation effects over the surface, in general. On impact, the geometrical changes in hoplengths from surface perturbations are given by (with reference to Fig. 2b)

$$ \Delta\ell(\ell, x_\downarrow) = \cot\theta(\ell)[Z(x_\uparrow) - Z(x_\downarrow)], \qquad (M5) $$

for impact angle $\theta(\ell)$. The behavior of geometrical $\Delta\ell$ over a modulated surface is shown in Figure M2. These geometrical effects on $\Delta\ell$ tend to cause impact focusing around the upwind inflection ripple point and is the main contributor to the flux modulations and $\ell \sim \lambda$ amplification—as illustrated in Figure 2b in the main text. The concept of hoplength changes $\Delta\ell$ has been used in past impact ripple models [7, 33], most notably in the limit of very small hoplengths where $\Delta\ell \sim -\cot\theta\ell\partial_x Z$ [7], but never in its general form as derived here. We also note that $\cot\theta(\ell)$ is another function measured in the DEM that shows approximately power-law like behaviour for the majority of particles in transport $\ell \lesssim \langle\ell\rangle$ (Fig. M2).

## Linking Ejections to Impacts

Given that Equation (M3) expresses the impact flux $\varphi_\downarrow$ in terms of the ejection flux $\varphi_\uparrow$ and hoplength modulations $\Delta\ell$ (Eq. M5), to close the problem we must relate ejections to impacts. We do this using a novel nonlocal splash-process approach. For an ensemble of impacts of a given hoplength $\ell$ at impact site $x_\downarrow$ there are on average $n_\epsilon(\epsilon|\ell)$ grains that are ejected at distance $\epsilon = x_\uparrow - x_\downarrow$ away from the impact site, per unit distance (as seen in the DEM (Fig. 2c) and experiments [34]). The rate of ejections and rebounds is given by

$$ \begin{aligned} \varphi_\uparrow(x_\uparrow) &= \iint n_\epsilon(\epsilon|\ell)\psi(\ell, x_\uparrow - \epsilon)d\epsilon d\ell \\ &\approx \int n_\uparrow(\ell)\psi(\ell, x_\uparrow)d\ell - \partial_{x_\uparrow} \int \bar{\epsilon}(\ell)n_\uparrow(\ell)\psi(\ell, x_\uparrow)d\ell, \end{aligned} \qquad (M6) $$

where the total number of ejections and rebounds per impact (the replacement capacity) is $n_\uparrow(\ell) = \int n_\epsilon(\epsilon|\ell)d\epsilon$ and the average impact-ejection lag distance per impact is $\bar{\epsilon}(\ell) = n_\uparrow(\ell)^{-1} \int \epsilon n_\epsilon(\epsilon|\ell)d\epsilon$. In general $n_\uparrow(\ell)$ and $\bar{\epsilon}(\ell)$ are functions of impact energy and other quantities, however, in the DEM we find the hoplength to be a good proxy for impact energy.
Equation (M6), evaluated in the homogeneous steady state, implies that the replacement capacity $n_↑(ℓ)$ is constrained to be

$$ \langle n_↑ \rangle = 1. \qquad (M7) $$

This useful constraint indicates that $n_↑(ℓ)P_0(ℓ) = P_n(ℓ)$ can be thought of as a probability distribution (SI II).

The constraint that $\langle n_↑ \rangle = \int n_↑(ℓ)P_0(ℓ)dℓ = 1$ for all steady-state transport conditions also implies that $\langle \varepsilon \rangle = \int \bar{\varepsilon}(ℓ)P_0(ℓ)dℓ$ varies little across transport conditions. Figure 3 suggests that $n_↑ \sim \bar{\varepsilon}^2$ for large density ratios $\rho_p/\rho_f$ (Mars, Earth) suggesting that $\bar{\varepsilon}^2$ is the (2D) volume of grains ejected from the surface. And for small $\rho_p/\rho_f$ (Venus), $n_↑ \sim \bar{\varepsilon}$ representing a more shear-like dislodging of surface grains. This behavior is characteristic of aeolian saltation and bedload transport, respectively [35]. Such relationships between $n_↑(ℓ)$ and $\bar{\varepsilon}(ℓ)$ suggests that can that the average $\langle \varepsilon \rangle \propto d$ is roughly constant across planetary conditions—which is what we find in the DEM (see caption of Fig. 3).

## Linear Instability: Wavelength and Speed Selection

Equations (M1-M6) provide a closed form model for surface evolution, given the flatbed quantities $P_0(ℓ)$, $n_↑(ℓ)$, $\bar{\varepsilon}(ℓ)$ and $\cot \theta(ℓ)$. To understand the dynamics of ripple emergence from a flat surface, we perform a linear stability analysis of the flat bed solution. Small perturbations to the surface $Z(x, t)$, such that $|\partial_x Z| \ll 1$, can be decomposed in Fourier space, $Z_k(x, t) = \hat{Z}(t)e^{ikx}$ (SI II B). The mode $k = 2\pi/\lambda$ grows exponentially in time according to

$$ \hat{Z}(t) \sim \hat{Z}(0)e^{(\sigma - ikc)t} \qquad (M8) $$

where $\sigma(k)$ is the growth rate and $c(k)$ the propagation speed.

The Fourier transform of the impact flux density obeys at the linear order

$$ \hat{\psi}(\ell, t) \sim \psi_0(\ell)(\mathcal{R} + i\mathcal{S}) k \hat{Z}(t) \qquad (M9) $$

for $\mathcal{R}(ℓ)$ and $\mathcal{S}(ℓ)$ are the in-phase and out-of-phase components the impact flux amplifications respectively, relative to the modulated surface, and are approximately unique functions of $ℓ/\lambda$—i.e., scale invariant (Fig.’s 2b and M3). The growth rate and speed of a given perturbation mode deduce as:

$$ \begin{aligned} \frac{\sigma}{k\varphi_0} &\sim \langle (1 - n_\uparrow)\mathcal{R} \rangle - k\langle \varepsilon_{\Delta\ell} \rangle_n \langle n_\uparrow\mathcal{S} \rangle \\ \frac{c}{\varphi_0} &\sim \langle (n_\uparrow - 1)\mathcal{S} \rangle \end{aligned} \qquad (M10) $$

where $\langle \varepsilon_{\Delta\ell} \rangle_n$ is a weighted average of $\bar{\varepsilon}(ℓ)$ over hoplength changes $\Delta\ell(ℓ)$—i.e., for geometrical hoplength changes (Eq. M5) $\langle \varepsilon_{\Delta\ell} \rangle_n = \langle \varepsilon_\theta \rangle_n = \langle \bar{\varepsilon} \cot \theta \rangle_n / \langle \cot \theta \rangle_n$ where $\langle \cdot \rangle_n = \int (\cdot) P_n dℓ$ (see SI Section II). This way of stating the
dispersion relation has a particularly intuitive meaning: The destabilizing mechanism represents the average number of grains that remain on the crest after traveling from the upwind face $\langle(1 - n_\uparrow)R\rangle$. This effect is countered by the stabilization mechanism from the increase number of grains being eroded from the crest due to the down wind shift in the ejection flux from $k\langle\epsilon_{\Delta\ell}\rangle_n\langle n_\uparrow S\rangle$. And the propagation speeds are determined by the number of grains eroded from or deposited on the upwind face $\langle(n_\uparrow - 1)S\rangle$ which is positive in the regime where the wavelength $\lambda$ is of order $\langle\ell\rangle$ because grains are being distributed from the upwind face over to the back of the ripple, and negative for wavelengths in the regime $\lambda \gg \langle\ell\rangle$.

The fastest growing wavelength of Equation (M10) is found at the maximum $\partial_k\sigma = 0$ which gives the scaling approximation provided in Equation (1) in the main text. Note that in Equation (1) in the main text we use $\langle\epsilon_{\Delta\ell}\rangle_n \equiv \langle\epsilon\rangle$ for simplicity. The wavelength and speed scalings from this model (Eq.'s M10 and 1) can be seen in Figure M4 over a wide range of transport conditions.

DEM Simulations

The discrete element model (DEM) used in this study to simulate sediment transport and ripple formation resolves grain-fluid and grain-grain forces on each individual particle. The simulations are quasi-3D with an along wind direction, $x$, of length 3000d with periodic boundaries, a vertical direction, $z$, of height 800d (higher than the largest hopheight), and a cross wind direction $y = 1d$. The domain is filled with ~30,000 grains of diameter $d$ with minor polydispersity. See Duran et al., 2012 [14], for a detailed description of the model dynamics.

In general, sediment transport is driven by fluid drag on sand grains. Thus the system is chiefly controlled by six parameters: the particle density $\rho_p$ and diameter $d$, the fluid density $\rho_f$, kinematic viscosity $\nu$, and shear velocity $u_*$, and gravity $g$. With these sediment transport can be characterized by three dimensionless numbers,

$$ \frac{\rho_p}{\rho_f}, \quad \text{Ga} = \frac{\sqrt{(\rho_p/\rho_f - 1)gd^3}}{\nu}, \quad \Theta = \frac{u_*^2}{(\rho_p/\rho_f - 1)gd}, \qquad (\text{M11}) $$

which are used to vary the transport conditions in the DEM. The density ratio between particles and the fluid $\rho_p/\rho_f$ is used to set the effective planetary regime, the Galileo number $Ga$ is used to set the grain size $d$, and the Shields number $\Theta$ sets the wind speed. Table M1 shows the values of these three numbers used in the DEM simulations presented in this study.

In our planetary simulations, ripples emerge spontaneously from flatbed, steady-state transport as seen in the space-time profiles in Figure 4b and in the Supplemental Video. Additionally, for each flatbed simulation we also conducted a range of simulations with the initial bed modulated as a sine wave of a given wavelength $\lambda$ and small amplitude ($\lesssim 3d$, much smaller than the imposed wavelength $\lambda$) to directly simulate the linear instability and measure the growth rates $\sigma(k)$ and speeds $c(k)$ for various perturbation modes $k = 2\pi/\lambda$ (as was done in [9]). As seen in Figure 4, the fastest growing wavelength measured by direct simulations of $\sigma(k)$ agrees with the wavelength that
spontaneously emerges from the flatbed simulations (diamond markers).

In order to obtain a closed analytical solution for $\sigma(k)$ and $c(k)$ based on our model derived above we approximate the flatbed distribution of hoplengths $P_0(\ell)$ as a power law with an exponential tail

$$P_0 \sim \ell^{-1+\alpha} e^{-\ell/\mu}, \qquad (M12)$$

where $\alpha$ takes into account the small deviation from the $\sim \ell^{-1}$ scaling observed for different density ratios (Fig. 3) and $\mu \sim \langle \ell \rangle$ sets the scale of the distribution—values of $\mu$ are chosen such that the analytical values of $\langle \ell \rangle$ matches the DEM measurements. For the normalization of $P_0$, $\int_{\ell_i}^{\infty} \ell^{-1+\alpha} e^{-\ell/\mu} d\ell$ we use a minimum hoplength of $\ell_i = 0.1d$ for $\alpha \le 0$ and $\ell_i = 0$ for $\alpha > 0$. We also approximate the other remaining functions used in the analytical model as power laws of the form

$$n_{\uparrow} \sim \ell^{\gamma_n}/\langle\ell^{\gamma_n}\rangle, \quad \bar{\varepsilon} \sim \bar{\varepsilon}_{1d}(\ell/d)^{\gamma_{\varepsilon}}, \quad \cot\theta \sim \cot\theta_{1d}(\ell/d)^{\gamma_{\theta}} \qquad (M13)$$

where $\langle \ell^{\gamma_n} \rangle$ comes from the constraint $\langle n_{\uparrow} \rangle = 1$. The values $\bar{\varepsilon}_{1d}$ and $\cot\theta_{1d}$ are the quantities when the hoplength $\ell = 1d$. Parameters used for calculation of the dispersion relation predictions in Figure 4 are seen in Table M1.

Direct measurements of $\bar{\varepsilon}(\ell)$ in the DEM during transport is complicated by the fact that we must measure $\bar{\varepsilon}$ for grains crossing a certain level (0.7d) above the average surface elevation and the fact the surface geometry is not smooth. The functional form of $\bar{\varepsilon}$ measurements seem to be robust when compared to single particle collision experiments (not shown) but measurement error is observed during transport, mostly for small hoplengths. Thus we fix $\bar{\varepsilon}_{1d} = 1d$ for all dispersion relation predictions—a physically reasonable constraint—and just fit the power $\gamma_{\varepsilon}$ (Table M1). The other parameters used in the model are approximated from the DEM measurements seen in Figures 3 and M2.

We lastly note that the model predictions of the dispersion relation ($\sigma$ and $c$) have been scaled by a proportionality factor to provide best fits for DEM measurements (Fig. 4). We do not expect the DEM measured dispersion relation to match the predictions exactly for three reasons: 1) the vertical flux scale $\phi_0$ is measured in the DEM at 0.7d above the surface and should therefore be smaller than the theoretical value of $\phi_0$ at the surface, 2) the analytical fitting functions of DEM quantities in Equations (M12) and (M13) are approximations, and 3) the hoplength changes $\Delta\ell$, which set the flux modulation magnitude (and therefore growth rate and speed magnitudes; SI II B), used in our analytical predictions consist of only geometrical contributions (Eq. M5) and are therefore smaller in magnitude than what is measured in the DEM. Namely, the DEM also simulates $\Delta\ell$ contributions from ejection angles responses to surface geometry—which would increase the overall magnitude of $\Delta\ell$ and add increasing flux modulation effects for small $\lambda$ (higher slopes). These additional effects are not the main contributors of the flux modulations (Fig. M3) and are thus effectively ignored, up to a constant, in this version of the analytical model. For these three reasons we introduce a multiplicative constant to fit the predicted dispersion relations for planetary condition simulated in the DEM. This constant is labeled as $\sigma_{\text{DEM}}/\sigma_{\text{pred}}$ in Table M1. Future work will be done to better understand
<table><caption>TABLE M1. DEM simulation parameters and associated fit values for the analytical model. The values of $Ga$ are chosen to keep $d \simeq 100\mu m$ and the values of $Θ$ are chosen such that the steady-state number of particles in transport per unit area is $\sim 1/50d^2$—representing fairly rarefied transport conditions.</caption><thead><tr><th></th><th>ρ<sub>p</sub>/ρ<sub>f</sub></th><th>Ga</th><th>Θ</th><th>μ/d</th><th>α</th><th>γ<sub>n</sub></th><th>γ<sub>ε</sub></th><th>γ<sub>θ</sub></th><th>ε&#773;<sub>1d</sub>/d</th><th>cot θ<sub>1d</sub></th><th>σ<sub>DEM</sub>/σ<sub>pred</sub></th></tr></thead><tbody><tr><td>Mars</td><td>10<sup>5</sup></td><td>2</td><td>0.017</td><td>4800</td><td>-0.15</td><td>0.52</td><td>0.27</td><td>0.24</td><td>1</td><td>0.95</td><td>0.29</td></tr><tr><td>Earth</td><td>2000</td><td>10</td><td>0.017</td><td>750</td><td>0</td><td>0.4</td><td>0.25</td><td>0.24</td><td>1</td><td>0.95</td><td>0.51</td></tr><tr><td>Titan</td><td>200</td><td>14</td><td>0.032</td><td>393</td><td>0</td><td>0.35</td><td>0.23</td><td>0.24</td><td>1</td><td>1.5</td><td>0.54</td></tr><tr><td>Venus</td><td>50</td><td>40</td><td>0.036</td><td>63</td><td>0.38</td><td>0.18</td><td>0.15</td><td>0.07</td><td>1</td><td>3</td><td>0.29</td></tr></tbody></table>

the consequences more precise analytical approximations (Eq.'s M12 and M13) and of additional contributions to hoplength modulations $\Delta\ell$ from ejection angle and wind modulation effects.

**Author Contributions:** All authors contributed to experimental design, interpretation and writing.

**Data Availability Statement:** Data will be available in an external repository after paper acceptance.

**Code Availability Statement:** Code for the discrete element model (DEM) used in this study can be made available upon request.

**Competing Interest Statement:** The authors declare no competing interests.
<figure><img src="image_5.png" /></figure>

FIG. M1. (a-d) DEM simulations of transport on Earth showing the time series of the number of grains in transport per unit area $n(t)$ rescaled by the average $N \sim 1/50d^2$. Panel (b) is the zoomed in box outlined in (a) and (c) is the black box in (b) and similarly (d) is the zoomed box of (c). Panel (e) shows the spectrum of fluctuations in $n(t)$ of a given frequency $|\hat{n}(f)|$ (i.e., the square root of the power spectrum). The spectrum shows power law decay as $1/f$ hinting at the critical sensitivity of the granular surface: With the wind speed at the bed buffered to the threshold of motion [9, 21], the bed is sensitive to perturbations. In this condition, when an impact ejects grains that then eject other grains, and so on, this feedback can cascade into large bursts of transport. These cascades cause fluctuations in transport across a range of scales. Highlighted by the solid markers in (e) are the frequency windows used in (a-d). Note that (b-d) fall in the scale-free region. The insets in (e) show the spectrum for Mars (red), Earth (blue), Titan (green) and Venus (gold) all showing similar scale-free $|\hat{n}| \sim 1/f$ spectra.
<figure><img src="image_6.png" /><figcaption>FIG. M2. (a) Cotangent of impact angles $cot θ(ℓ)$ for Mars (red), Earth (blue), Titan (green) and Venus (gold) conditions, measured from DEM simulations. $cot θ(ℓ)$ is rescaled by its average for clarity of presentation. For Mars, Earth, Titan and Venus $⟨cot θ⟩ \simeq (1.4, 1.5, 2.8, 3.7)$, respectively. (b) The Fourier transform of geometric hoplength modulations (for small surface perturbations $\hat{Z}$) are defined as $\Delta\hat{\ell} = cot θ(ℓ)(e^{-ikℓ} - 1)\hat{Z}$ (Eq. M5). Plotted in (b) is the (negative) phase and magnitude of $\Delta\hat{\ell}/\hat{Z}$ versus rescaled hoplength $ℓ/λ$ for $λ = 200d$ in Earth conditions ($⟨ℓ⟩ \simeq 100d$). Also shown over the magnitude is the envelop $2cot θ(ℓ)$.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>FIG. M3. Hoplength amplification rates $R(kℓ)$ and $S(kℓ)$ rescaled by $⟨cot θ⟩_n = ∫ cot θ(ℓ)P_n(ℓ)dℓ$ (SI II F). Circles are DEM measurements for Earth conditions and lines are predictions from the analytical model using geometrical hoplength changes (Fig. M2). DEM markers are scaled down by a factor of 1.5 to take into account the missing $\Deltaℓ$ contributions. The dashed lines are the approximations $R \propto \sin kℓ$ and $S \propto 1 + \cos kℓ$ (see SI II E).</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>FIG. M4. Dispersion relation solutions over a wide range of planetary transport conditions (represented by changing $\langle l \rangle$). Positive values are represented by solid lines and negative values by dashed lines. The rescaled growth rate $\sigma$ versus $\lambda$ rescaled by the weighted average lag length $\langle \epsilon_{\theta} \rangle_n = \langle \epsilon \cot \theta \rangle_n / \langle \cot \theta \rangle_n$ shows that the fastest growing wavelength in this model scales with $\langle \epsilon_{\theta} \rangle_n$. The rescaled speed $c$ versus $\lambda$ rescaled by average hoplength $\langle l \rangle$ shows that the ripple propagation speed in this model is a unique function of $\lambda/\langle l \rangle$—and when evaluated at the fastest growing wavelength $\langle \epsilon_{\theta} \rangle_n \propto d$ the associated speed is $c/\varphi_0 = f(\langle l \rangle/d)$. These solutions were derived for $\alpha = 0, \gamma_n = 0.4, \gamma_{\epsilon} = 0.25, \gamma_{\theta} = 0.25$ and $\tilde{\epsilon}_{1d}/d = \cot \theta_{1d} = 1$. These parameters represent typical conditions for saltation, $\rho_p/\rho_f \gtrsim 100$ (Table M1). We can see that for $\lambda \sim 200d$ the ripple-antiripple transition ($c=0$) occurs at a critical hoplength $\langle l \rangle_c \sim 60d$ consistent with stationary ripples simulations observed for Titian (Fig. 4). To a good approximation the speeds evaluated at the fastest growing wavelength scale as $c/\varphi_0 \propto (\langle l \rangle/\langle l \rangle_c)^{2/5} - 1$ for the parameters used here.</figcaption></figure>
[1] F. Charru, B. Andreotti, and P. Claudin, Sand ripples and dunes, Annual Review of Fluid Mechanics 45, 469 (2013).

[2] O. Duran Vinent, B. Andreotti, P. Claudin, and C. Winter, A unified model of ripples and dunes in water and planetary environments, Nature Geoscience 12, 345 (2019).

[3] R. A. Bagnold, The physics of blown sand and desert dunes (Courier Corporation, 2012).

[4] M. Lapotre, R. Ewing, M. Lamb, W. Fischer, J. Grotzinger, D. Rubin, K. Lewis, M. Ballard, M. Day, S. Gupta, et al., Large wind ripples on mars: A record of atmospheric evolution, Science 353, 55 (2016).

[5] H. Yizhaq, K. Tholen, L. Saban, N. Swet, C. Lester, S. Silvestro, K. R. Rasmussen, J. P. Merrison, J. J. Iversen, G. Franzese, et al., Coevolving aerodynamic and impact ripples on earth, Nature Geoscience, 1 (2024).

[6] R. S. Anderson, A theoretical model for aeolian impact ripples, Sedimentology 34, 943 (1987).

[7] L. Prigozhin, Nonlinear dynamics of aeolian sand ripples, Physical Review E 60, 729 (1999).

[8] Z. Csahók, C. Misbah, F. Rioual, and A. Valance, Dynamics of aeolian sand ripples, The European Physical Journal E 3, 71 (2000).

[9] O. Durán, P. Claudin, and B. Andreotti, Direct numerical simulations of aeolian sand ripples, Proceedings of the National Academy of Sciences 111, 15665 (2014).

[10] B. Andreotti, P. Claudin, J. J. Iversen, J. P. Merrison, and K. R. Rasmussen, A lower-than-expected saltation threshold at martian pressure and below, Proceedings of the National Academy of Sciences 118, e2012386118 (2021).

[11] J. Miller, J. Marshall, R. Greeley, et al., Wind ripples in low density atmospheres, NASA Technical Memorandum 89810, 268 (1987).

[12] B. Andreotti, P. Claudin, and O. Pouliquen, Aeolian sand ripples: Experimental study of fully developed states, Physical review letters 96, 028001 (2006).

[13] R. Greeley, J. Marshall, and R. Leach, Microdunes and other aeolian bedforms on venus: Wind tunnel simulations, Icarus 60, 152 (1984).

[14] O. Durán, B. Andreotti, and P. Claudin, Numerical simulation of turbulent sediment transport, from bed load to saltation, Physics of Fluids 24 (2012).

[15] T. Pähtz and O. Durán, Unification of aeolian and fluvial sediment transport rate from granular physics, Physical Review Letters 124, 168001 (2020).

[16] T. Pähtz and O. Durán, Scaling laws for planetary sediment transport from dem-rans numerical simulations, Journal of Fluid Mechanics 963, A20 (2023).

[17] T. D. Ho, A. Valance, P. Dupont, and A. O. El Moctar, Aeolian sand transport: Length and height distributions of saltation trajectories, Aeolian Research 12, 65 (2014).

[18] K. R. Rasmussen, A. Valance, and J. Merrison, Laboratory studies of aeolian sediment transport processes on planetary surfaces, Geomorphology 244, 74 (2015).

[19] M. Lämmel and K. Kroy, Analytical mesoscale modeling of aeolian sand transport, Physical Review E 96, 052906 (2017).

[20] G. A. Held, D. Solina, H. Solina, D. Keane, W. Haag, P. Horn, and G. Grinstein, Experimental study of critical-mass fluctuations in an evolving sandpile, Physical Review Letters 65, 1120 (1990).

[21] O. Durán, P. Claudin, and B. Andreotti, On aeolian transport: Grain-scale interactions, dynamical mechanisms and scaling
laws, Aeolian Research 3, 243 (2011).

[22] L. Yan, A. Barizien, and M. Wyart, Model for the erosion onset of a granular bed sheared by a viscous fluid, Physical Review E 93, 012903 (2016).

[23] E. T. Owens and K. E. Daniels, Sound propagation and force chains in granular materials, Europhysics Letters 94, 54005 (2011).

[24] J. D. Walker, An experimental study of wind ripples, Ph.D. thesis, Massachusetts Institute of Technology (1981).

[25] E. Schmerler, I. Katra, J. F. Kok, H. Tsoar, and H. Yizhaq, Experimental and numerical study of sharp's shadow zone hypothesis on sand ripple wavelength, Aeolian research 22, 37 (2016).

[26] H. Cheng, C. Liu, J. Li, B. Liu, Z. Zheng, X. Zou, L. Kang, and Y. Fang, Experimental study of aeolian sand ripples in a wind tunnel, Earth Surface Processes and Landforms 43, 312 (2018).

[27] M. A. Makeev, R. Cuerno, and A.-L. Barabasi, Morphology of ion-sputtered surfaces, Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms 197, 185 (2002).

[28] S. A. Norris and M. J. Aziz, Ion-induced nanopatterning of silicon: Toward a predictive model, Applied Physics Reviews 6, 10.1063/1.5043438 (2019).

[29] B. Werner and D. Gillespie, Fundamentally discrete stochastic model for wind ripple dynamics, Physical review letters 71, 3230 (1993).

[30] W. Landry and B. Werner, Computer simulations of self-organized wind ripple patterns, Physica D: Nonlinear Phenomena 77, 238 (1994).

[31] D. A. Huntley, G. Coco, K. R. Bryan, and A. B. Murray, Influence of “defects” on sorted bedform dynamics, Geophysical Research Letters 35 (2008).

[32] B. Werner and G. Kocurek, Bedform spacing from defect dynamics, Geology 27, 727 (1999).

[33] H. Nishimori and N. Ouchi, Formation of ripple patterns and dunes by wind-blown sand, Physical Review Letters 71, 197 (1993).

[34] M. Ammi, L. Oger, D. Beladjine, and A. Valance, Three-dimensional analysis of the collision process of a bead on a granular packing, Physical Review E 79, 021305 (2009).

[35] T. Pähtz and O. Durán, Fluid forces or impacts: What governs the entrainment of soil particles in sediment transport mediated by a newtonian fluid?, Physical Review Fluids 2, 074303 (2017).

[36] J. K. Pierce and M. A. Hassan, Joint stochastic bedload transport and bed elevation model: Variance regulation and power law rests, Journal of Geophysical Research: Earth Surface 125, e2019JF005259 (2020).

[37] P. Sigmund, A mechanism of surface micro-roughening by ion bombardment, Journal of Materials Science 8, 1545 (1973).

[38] D. J. Furbish, P. K. Haff, J. C. Roseberry, and M. W. Schmeeckle, A probabilistic description of the bed load sediment flux: 1. theory, Journal of Geophysical Research: Earth Surface 117, 10.1029/2012JF002352 (2012).
# Supplemental Information: Emergence of wind ripples controlled by mechanics of grain-bed impacts

## I. SURFACE EVOLUTION

### A. The Surface

Our model is based on the fundamental constraint that sand grains are conserved. We thus attempt to evolve the surface $Z(x, t)$ based on this mass conservation constraint, focusing on one spatial dimension $x$ in the direction of the wind (generalization to two spatial dimensions $(x, y)$ is straight forward). However, the surface elevation $Z(x, t)$ is not trivially defined [36]. In the simplest case, the surface could be defined by the elevation below which the volume fraction of grains is constant, at the bed packing fraction $\phi_b$. We will call this elevation $Z_b(x, t)$. In practice $\phi_b$ and $Z_b$ are averaged over a finite surface area $dxdy$ (here $dy = 1d$). For example, if we want to observe elevation changes on the scale of a single grain of diameter $d$, such that $dxdy = d^2$, then the surface slopes would change drastically and discontinuously over space. Thus we seek to coarse-grain the surface over an area much larger than a single grain but not too large to where we are missing important changes in surface slope. In practice we average the surface over $10d$ in DEM measurements.

The elevation $Z_b(x, t)$ becomes less identifiable as the surface as the transport rate $Q$ increase. This is because the volume fraction of grains at a given elevation $\phi(z)$ is no longer a simple step-function between the bed $\phi_b$ and above due to bed dilation from grain-bed impacts. Thus we consider the evolution of the effective surface given transport, $Z(x, t)$, which is defined as

$$Z(x, t) = \int_{0}^{\infty} \frac{\phi(z, x, t)}{\phi_b} dz. \qquad (S1)$$

If $M(x, t)$ is the total volume of grains in motion per unit area—proportional to the shear stress $\tau$ above the transport threshold $M \sim \tau - \tau_{th}$ [15]—then

$$Z(x, t) = Z_b(x, t) + \frac{M(x, t)}{\phi_b}. \qquad (S2)$$

Note that in steady state transport over a flat bed, $Z(x)$ must on average be equal to the elevation of the surface without transport $Z_b(M = 0)$ because $Z$ simply represents the elevation if all grains in transport were piled on the surface $Z_b$. Thus, $Z_b(\tau)$ must decrease at the same rate as $M(\tau)$ increases with $\tau - \tau_{th}$.
## B. Mass Conservation

As we have to coarse grain the spatial scale to define the surface elevation $Z(x, t)$, we also coarsen the time scale for surface evolution $dt$ to be slow compared to granular time scales—for example, $\sqrt{d/g}$ being (roughly) the time it takes for a grain to fall a distance $d$ under gravity $g$. I.e., if $dt = \sqrt{d/g}$ then the evolution of surface height over time $\partial_t Z(x, t)$ would change discontinuously by large amounts. So we take $dt$ to be much larger than $\sqrt{d/g}$ but much smaller than the ripple growth time scale $\sigma^{-1}$, defined by the growth rate $\sigma$ which scales with the excess shear stress $\tau - \tau_{th}$. For reference, the DEM suggests that for Earth-like conditions $\sigma^{-1} \sim \mathcal{O}(10^5 \sqrt{d/g})$ for rarefied transport (e.g., for one particle in transport per $50d^2$) and the average travel time for a grain trajectory of hoplength $\ell \sim 100d$ is $\mathcal{O}(10\sqrt{d/g})$.

Given the above considerations, the equation of motion for the surface $Z(x, t)$ is given by the gradient in the total horizontal flux $Q(x, t)$ (normalized by $\phi_b$) defined as the volume of grains crossing a vertical plain per unit width and time:

$$ \partial_t Z(x, t) = -\partial_x Q(x, t). \qquad (S3) $$

This is the so-called “Exner Equation”. Whereas the evolution of the static bed $Z_b(x, t)$ is given by the difference in the vertical fluxes of grains impacting the surface and leaving the surface, $\varphi_\downarrow(x, t)$ and $\varphi_\uparrow(x, t)$ (again, all normalized by $\phi_b$):

$$ \partial_t Z_b(x, t) = \varphi_\downarrow(x, t) - \varphi_\uparrow(x, t). \qquad (S4) $$

We assume the change in the transport “activity” above the static surface over time to be negligible $\partial_t M \simeq 0$ in comparison to $\partial_t Z_b$. This is an argument that the time it takes for transport to adjust is much faster than the exceedingly slow time scale evolution of the surface. Thus surface evolution is given by

$$ \partial_t Z(x, t) = -\partial_x Q(x, t) = \varphi_\downarrow(x, t) - \varphi_\uparrow(x, t). \qquad (S5) $$

## C. Linking Impacts to Ejections

Our model for ripple growth is derived by relating fluxes at each point to upwind fluxes using grain trajectories. The impact flux $\varphi_\downarrow(x_\downarrow)$ is found by relating all grains impacting a location $x_\downarrow$ to the location they left from $x_\uparrow$—linking $\varphi_\downarrow(x_\downarrow)$ to $\varphi_\uparrow(x_\uparrow)$ via the hoplength $\ell = x_\downarrow - x_\uparrow$. This requires information about the distribution of hoplengths, which is generally modulated in response to a modulated surface $Z(x, t)$.

We define $\psi(\ell, x_\downarrow)d\ell dx_\downarrow$ as the volume of grains (at bed packing fraction $\phi_b$) impacting the surface between $x_\downarrow$ and $x_\downarrow + dx_\downarrow$ after hoplengths between $\ell$ and $\ell + d\ell$, such that $\varphi_\downarrow(x_\downarrow) \equiv \int \psi(\ell, x_\downarrow)d\ell$. The conservation of mass, with
reference to Figure 2a,b in the main text, dictates

$$ \psi(\ell, x_{\downarrow})d\ell dx_{\downarrow} = \varphi_{\uparrow}(x_{\uparrow})dx_{\uparrow}P_{0}(\ell_{0})d\ell_{0}, \qquad (S6) $$

which states that the grains impacting at $x_{\downarrow}$ with modulated hoplengths $\ell$ are the same grains that left from $x_{\uparrow}$ with unmodulated hoplengths $\ell_0$ given by the homogeneous probability density $P_0(\ell_0)$ defined on a flat surface.

In order to equate the impact flux density $\psi(\ell, x_{\downarrow})$ nonlocally to the upstream vertical flux $\varphi_{\uparrow}(x_{\uparrow})$ given the hoplength distribution $P_0(\ell_0)$ in Equation (S6), we change coordinates from $(\ell_0, x_{\uparrow})$-space associated with flat bed information to $(\ell, x_{\downarrow})$-space which takes into account the changes in the flat bed coordinates due to the modulated surface:

$$ \psi(\ell, x_{\downarrow}) = \varphi_{\uparrow}(x_{\downarrow} - \ell)P_{0}(\ell - \Delta\ell)|\mathcal{J}_{\ell, x_{\downarrow}}|, \qquad (S7) $$

where $\Delta\ell(\ell, x_{\downarrow}) = \ell - \ell_0$ is the small perturbation to the flat bed hoplength $\ell$ from surface modulation effects and $|\mathcal{J}_{\ell, x_{\downarrow}}| = |1 - \partial_{\ell}\Delta\ell - \partial_{x_{\downarrow}}\Delta\ell|$ is the resulting Jacobian of the coordinate transformation.

## D. Linking Ejections to Impacts

Equation (S7) puts the impact flux $\varphi_{\downarrow}$ in terms of the ejection flux $\varphi_{\uparrow}$ nonlocally using grain trajectories $\ell$. To close the problem we now relate ejections to impacts—which is a more local response. Reiterating the ejection process outlined in the Methods, ejections and rebounds are generally described using the replacement capacity density (per unit lag distance $\varepsilon = x_{\uparrow} - x_{\downarrow}$), $n_{\varepsilon}(\varepsilon|\ell, x_{\downarrow})$—which may in general depend on surface geometry and thus may be a function of space. Thus the general rate of ejections and rebounds is given by

$$ \varphi_{\uparrow}(x_{\uparrow}) = \iint \psi(\ell, x_{\uparrow} - \varepsilon)n_{\varepsilon}(\varepsilon|\ell, x_{\uparrow} - \varepsilon)d\varepsilon d\ell. \qquad (S8) $$

Where the replacement capacity (total number of ejections/rebounds per impact) is expressed as $n_{\uparrow}(\ell, x) = \int n_{\varepsilon}(\varepsilon|\ell, x)d\varepsilon$. A spatial dependence of $n_{\varepsilon}$ suggests that there are two modulating effects on the fluxes, one from modulations of trajectories $\Delta\ell$ and the other from modulations of ejections and rebounds caused by impact-ejection interactions with surface geometry. The latter is the main driver for the modulation of surface erosion and the emergence of ion sputter ripples [28, 37] where ion impacts penetrate into a material and sputter (erode) the surface in a way that greatly depends on surface slope and curvature. Although this is a mechanism similar to the one in which grains are ejected from a sand bed, we find that this effect is in fact small compared to the modulation component coming from $\Delta\ell$—which is reasonable given that grain impacts do not travel deep into the bed surface in comparison to the ion-sputter process [37]. Thus from here on we let the replacement capacity (density) be spatially
constant in our analysis,

$$n_{\varepsilon}(\varepsilon|\ell, x_{\downarrow}) \sim n_{\varepsilon}(\varepsilon|\ell). \qquad (S9)$$

A more general derivation of the model with spatially explicit ejection dependence is straight forward.

## II. LINEAR INSTABILITY

### A. Linear Regime

In this study we are mainly concerned with the linear instability of ripple formation from the uniform flat bed fixed point $Z_*(x) \equiv 0$. In the linear regime we know that perturbations to the surface $Z(x, t)$ must evolve as $\partial_t Z \sim \Omega Z$ where the growth rate $\Omega$ is a linear operator that determines how perturbation amplitudes grow vertically and propagate in space—i.e., $\Omega$ consists of local operations $\partial_x$, $\partial_{xx}$, $\partial_{xxx}$... and much more complicated spatial dependencies because the problem is generally nonlocal. The criterion for the linear regime is $|\partial_x Z| \ll 1$ for all $x$—for an average surface elevation over $dx$ bigger than a grain diameter $d$. Linearity allows us to explicitly expand space and time dependent quantities into a zeroth order component associated with flat-bed/steady-state transport and a first order component that takes into account small changes in proportion with $Z(x, t)$. Thus for example we let $\varphi(x, t) = \varphi_0 + \varphi_1(x, t)$ and $\psi(\ell, x, t) = \psi_0(\ell) + \psi_1(\ell, x, t)$. This allows the growth rate operator to be written explicitly as

$$\Omega = \frac{\varphi_{\downarrow 1}(x, t) - \varphi_{\uparrow 1}(x, t)}{Z(x, t)} \qquad (S10)$$

where we have used $\varphi_{\uparrow 0} = \varphi_{\downarrow 0} = \varphi_0$, which defines steady state on a flatbed. This shows that only first order perturbations to the fluxes are responsible for ripple growth. Because of this we will drop the subscript's-1 for all first order (spatially dependent) expressions below and explicitly denote flat bed (spatially constant) quantities with subscript's-0.

Imposing linearity, the first order perturbation to the impact flux density $\psi(\ell, x)$ in Equation (S7) simplifies to

$$\psi(\ell, x) = \varphi_{\uparrow}(x - \ell)P_0(\ell) - \varphi_0 \partial_{\ell, x} \left( \Delta\ell(\ell, x) P_0(\ell) \right) \qquad (S11)$$

where $\partial_{\ell, x} = \partial_{\ell} + \partial_x$. This is the same as Methods Equation (M2) and integrates to give the total impact flux $\varphi_{\downarrow}$ in Methods Equation (M3).

With Equation (S11) (and Methods Eq. M3), the surface evolution Equation (S5) can now be written to first order as

$$\partial_t Z(x, t) = \int \varphi_{\uparrow}(x - \ell, t) P_0(\ell) d\ell - \varphi_{\uparrow}(x, t) - \varphi_0 \partial_x \langle \Delta\ell(x, t) \rangle. \qquad (S12)$$
We see that the first term assumes no explicit changes to grain hoplengths—as has been commonly used in past impact ripple research [6, 8]—and the final term is a result of the extra flux associated with small changes in grain trajectories in the presence of the modulated surface.

Alternatively, using the definition of $\varphi_{\uparrow}$ from Equation (S8) the evolution of the surface may be described in terms of impacts and resulting bed response from ejections:

$$ \partial_t Z(x, t) = \iint \psi(\ell, x - \varepsilon, t) [\delta(\varepsilon) - n_{\varepsilon}(\varepsilon|\ell)] d\varepsilon d\ell, \qquad (S13) $$

where the net number of grains exchanged with the bed $\delta(\varepsilon) - n_{\varepsilon}(\varepsilon|\ell)$ is analogous to the “crater function” used in describing the growth of ion-sputter ripples [28].

## B. Growth and Modulation Rates

We can greatly simplify the analytical expressions above by independently analyzing perturbations of a given (Fourier) mode $k = 2\pi/\lambda$,

$$ Z_k(x, t) = \hat{Z}(k, t)e^{ikx} \qquad (S14) $$

where the sum over all $k$ recovers the surface $Z(x, t) = (2\pi)^{-1} \int Z_k(x, t)dk$, keeping the real part only. In the linear growth regime, $k|\hat{Z}| \ll 1$, the complex amplitude evolves according to $\hat{Z}(t) \propto e^{\Omega t}$ where the operator $\Omega(k)$ at a given mode $k$ turns into a complex growth rate,

$$ \Omega(k) = \sigma(k) + i\omega(k), $$

with $\sigma$ being the real growth rate of the magnitude $|\hat{Z}|$ and $c = -\omega/k$ is the propagation speed of the perturbed mode $k$.

At mode $k$ all other spatially explicit (first order) quantities can be cast as Fourier amplitudes proportional to $\hat{Z}$. We thus let

$$ \hat{\varphi}_{\downarrow} = \varphi_0 \mathcal{F}_{\downarrow} k \hat{Z} \quad \text{and} \quad \hat{\varphi}_{\uparrow} = \varphi_0 \mathcal{F}_{\uparrow} k \hat{Z} \qquad (S15) $$

which define the vertical flux modulation rates $\mathcal{F}_{\uparrow(\downarrow)}$ as complex numbers whose magnitude gives the size of the flux perturbation and whose phase locates the modulated flux relative to the surface phase. For example, $\arg \mathcal{F}_{\downarrow} = \pi/2$ says that the perturbation to the impact flux $\hat{\varphi}_{\downarrow}$ is located at the upwind inflection point—i.e., $\pi/2$ out of phase with
$\hat{Z}$. For the impact flux density and hoplength modulations we let

$$ \hat{\psi}(\ell) = \psi_0(\ell)\mathcal{A}(\ell)k\hat{Z} \quad \text{and} \quad \Delta\hat{\ell}(\ell) = \mathcal{L}(\ell)\hat{Z} \qquad (S16) $$

which defines the complex modulation rates $\mathcal{A}(\ell)$ and $\mathcal{L}(\ell)$, which encode the same information as $\mathcal{F}_{\uparrow(\downarrow)}$ but on a per hoplength basis. For example, $\mathcal{A}(\lambda)$ gives the amplification of the impact flux for grains with hoplength $\ell = \lambda$. We note that for geometrical hoplength changes (see Methods Eq. M5), the hoplength modulation rate is given by

$$ \mathcal{L}(\ell) = \cot\theta(\ell)(e^{-ik\ell} - 1), \qquad (S17) $$

where we used the identity $Z(x_\uparrow) = Z(x_\downarrow - \ell) = \hat{Z}e^{ikx_\downarrow}e^{-ik\ell}$.

Given these definitions we can simplify the surface evolution Equations (S5) and (S12) to find the so-called “dispersion relation” for the growth rate $\Omega(k)$:

$$ \frac{\Omega}{k\varphi_0} = \mathcal{F}_{\downarrow} - \mathcal{F}_{\uparrow} = (\tilde{P}_0 - 1)\mathcal{F}_{\uparrow} - i\langle\mathcal{L}\rangle \qquad (S18) $$

where $\tilde{P}_0 = \int P_0(\ell)e^{-ik\ell}d\ell$. The ascending flux modulation rate is defined by Equation (S8) as

$$ \mathcal{F}_{\uparrow} = \int \tilde{n}_{\varepsilon}(\ell)\mathcal{A}(\ell)P_{0}(\ell)d\ell, \qquad (S19) $$

where $\tilde{n}_{\varepsilon}(\ell) = \int n_{\varepsilon}(\varepsilon|\ell)e^{-ik\varepsilon}d\varepsilon$. Thus we are also able to write the dispersion relation in terms of the amplification of impacting trajectories and the resulting bed response (Eq. S13):

$$ \frac{\Omega}{k\varphi_0} = \int (1 - \tilde{n}_{\varepsilon}(\ell))\mathcal{A}(\ell)P_0(\ell)d\ell. \qquad (S20) $$

From the derivation of $\psi(\ell, x)$ in Equation (S11), the hoplength amplification rate is

$$ \mathcal{A} = \mathcal{F}_{\uparrow}e^{-ik\ell} - k^{-1}\{\mathcal{L}' + (P'_{0}/P_{0} + ik)\mathcal{L}\}, \qquad (S21) $$

where $(\cdot)' = \partial_\ell(\cdot)$. Thus we obtain a closed-form solution to $\Omega(k)$ by solving for $\mathcal{F}_{\uparrow}$ using Equations (S19) and (S21) to obtain

$$ \rho_n \mathcal{F}_{\uparrow} = -i \int \tilde{n}_{\varepsilon} \mathcal{L} P_0 d\ell + k^{-1} \int \tilde{n}_{\varepsilon}' \mathcal{L} P_0 d\ell, \qquad (S22) $$

where $\rho_n = 1 - \int \tilde{n}_{\varepsilon} P_0 e^{-ik\ell} d\ell$ is a factor associated with the condition that the vertical flux is in modulated equilibrium (Fig. 2b).
### C. Local and Nonlocal Regimes

It is useful to analyze the linear instability dynamics in terms of ‘local’ and ‘nonlocal’ flux response, which arise from the relationship between the perturbation size $\lambda$ and the average hoplength $\langle l \rangle$. In this system these regimes are identified via the function $1 - \tilde{P}_0$, which represents the modulated relationship between the flux modulation $\hat{Q}_0$ in response to a sinusoidal modulation of the ejection flux $\hat{\phi}_{\uparrow 0}$ on a flatbed (i.e., no $\Delta l$ contributions):

$$ \frac{1 - \tilde{P}_0}{ik} \equiv \frac{\hat{Q}_0}{\hat{\phi}_{\uparrow 0}} \sim \begin{cases} \langle l \rangle e^{-\gamma ik \langle l \rangle}, & \lambda/\langle l \rangle \gg 1 \\ \lambda e^{-i\pi/2}, & \lambda/\langle l \rangle \lesssim 1. \end{cases} \qquad (S23) $$

In the limit of long wavelength modulations relative to $\langle l \rangle$, the flux response is local and simply reduces to the identity $\hat{Q}_0 \sim \langle l \rangle \hat{\phi}_{\uparrow 0}$ with a downwind lag of $\hat{Q}_0$ relative to $\hat{\phi}_{\uparrow 0}$ proportional to $\langle l \rangle$ (Fig. S1)—with $\gamma$ being the proportionality constant in Eq. (S23). For wavelengths commensurate with $\langle l \rangle$, the phase relation between $\hat{Q}_0$ relative to $\hat{\phi}_{\uparrow 0}$ grows because of the nonlocal contributions from the fluxes upwind and asymptotically approaches a relative phase shift of $\pi/2$. This asymptotic phase shift behavior is purely due to the power-law-like shape of $P_0(l)$ which provides exceedingly large local contributions ($l/\lambda \lesssim 1$) to $\hat{Q}_0$ that dominate in the long-wavelength regime, but are counteracted by the increasing contributions from trajectories originating on other parts of a ripple ($l/\lambda \gtrsim 1$) as wavelength decreases relative to $\langle l \rangle$. And in the nonlocal limit the flux ratio magnitude $|\hat{Q}_0/\hat{\phi}_{\uparrow 0}| \sim \lambda$ because the nonlocal flux contributions from upwind scale with $\lambda$, almost irrespective of $P_0(l)$ (for small enough $\lambda/\langle l \rangle$).

### D. Salton-Repton Model

The seminal salton-repton approach to understanding wind ripple formation was first formalized by Anderson [6] and then later iterations of this model were proposed, all with similar concepts and predictions [7–9]—notably the prediction that ripple wavelengths are selected by a transport length scale, the average hoplength or repton length. In this framework, a uniform salton rain bombards the surface and causes ejections of reptons. Because the salton rain is uniform (unmodulated) and the impact-ejection lag length $\epsilon$ is ignored, the maximum in the ejection flux of reptons is located at the upwind inflection points where slopes are greatest. In our framework, this translates to $\phi_{\uparrow} \propto \partial_x Z$ or $F_{\uparrow} \propto i$ (Eq. S15) and thus

$$ \Omega_{\text{salt-rep}} \propto ik(\tilde{P}_0 - 1) \propto k^2 \frac{\hat{Q}_0}{\hat{\phi}_{\uparrow 0}} \qquad (S24) $$

where $P_0(l)$ can now effectively be interpreted and the distribution of repton hoplengths. In the salton-repton model the horizontal flux $Q(x)$ is shifted downwind of the ejection flux $\phi_{\uparrow}(x)$ by a distance depending only on the homogeneous probability distribution of hoplengths $P_0(l)$ (Fig. S2a,b)—and doesn’t take into account effects from
<figure><img src="image_9.png" /><figcaption>FIG. S1. Magnitude and phase relationships for the modulated flux $\hat{Q}_0$ responding to a sinusoidal ejection flux $\hat{\phi}_{\uparrow 0}$ at mode $k = 2\pi/\lambda$ on a flatbed. Markers are derived from DEM measurements of $P_0(\ell)$ for Mars, Earth, Titan and Venus conditions (see Fig. 3) and the solid lines are using power-law approximations of $P_0(\ell)$ with exponential cutoffs (Methods; Table M1).</figcaption></figure>

changes in hoplengths $\Delta\ell$ or the impact-ejection lag $\varepsilon$ (Fig. 2). And because of the power-law like structure of $P_0$ (Fig. 3), the salton-repton model predicts that the growth rate $\sigma_{\text{salt-rep}}(k) > 0$ for all $k$ because the phase shift $\arg\{\hat{Q}_0/\hat{\phi}_{\uparrow 0}\}$ doesn't allow for the flux maximum to pass the crest for $\hat{\phi}_{\uparrow} = \hat{\phi}_{\uparrow 0}$ maximized at the upwind infection point, $\pi/2$ out of phase (Fig. S1). In other words, the salton-repton model precludes erosion of the ripple crest for all $k$ (see Figures 4, S1 and S2).

It is often assumed that $P_0(\ell)$ is Gaussian-like, which is unphysical because ejections arise from a broad continuum of impacting grain energies (Fig. 3). Another common approach is to assume a single (reptating) grain trajectory $\ell = \langle\ell\rangle$ such that $P_0 \sim \delta(\ell - \langle\ell\rangle)$. Figure S2 shows the consequences of these unrealistic Gaussian-like and delta hoplength distributions when compared to a more realistic power-law like $P_0(\ell)$. Figure S2b shows the transient effects of a sudden step in the flatbed ejection flux $\phi_0$ at $x=0$ on the resulting relaxation of the horizontal flux $Q_0(x)$. Here the flatbed horizontal flux $Q_0(x)$ is the volume of grains ejected at $x-\ell$ that exceed the travel distance $\ell$ [38]. Thus $Q_0(x) = \phi_0 \int_0^x R(\ell)d\ell$ relaxes to its saturated state $\langle\ell\rangle\phi_0$ depending on the shape of $R(\ell) = \int^\ell P_0 d\ell$. For Gaussian/delta-like distributions this relaxation occurs quickly compared to $\langle\ell\rangle$ whereas for the power-law like
<figure><img src="image_10.png" /></figure>

FIG. S2. (a) Schematic of the widely used salton-repton model for ripple formation. A uniform salton rain is focused on upwind faces and causes ejections of reptons. Thus the maximum in the ejection flux $\varphi_{\uparrow}(x)$ (blue) is located at the upwind inflection points and the horizontal flux $Q(x)$ (red) is shifted downwind. Note that for the schematic drawn in (a), the flux $Q(x)$ is converging ($\partial_x Q < 0$) at the crest causing this particular ripple wavelength to grow. Section (b) on the right shows a set of grain trajectories being transported from the upwind inflection point. The thickness of each trajectory is given by the number of grains, or probability $P_0(\ell)$ with that given hoplength—taken to be Gaussian-like. The $P_0$ plot shows Gaussian-like (solid line), power-law like (dashed line) and delta (thin vertical line) distributions. The $Q_0(x)$ plot shows transient relaxation effects for a sudden step in the ejection flux to $\varphi_0$ at $x = 0$ (bottom plot, blue line) for each of the three distributions above it. Panel (c) shows the growth rates versus rescaled wavenumber $(\ell/\lambda)$ for the Gaussian-like $P_0$ (solid line; “unrealistic”) and power-law like $P_0$ (dashed line; “realistic”). We are also showing in the insets of (c) that the ratio $(\ell/\lambda)$ separates the dynamics from local $(\ell/\lambda \ll 1)$, where the form of $P_0$ has little importance, to nonlocal $(\ell/\lambda \gtrsim 1)$ where the form of $P_0$ is crucial. Red lines in the insets represent the flux $Q(x)$.

distribution the relaxation is much slower—saturating at $x \gg \langle\ell\rangle$. This means that the effects of perturbations to the vertical flux $\varphi_{\uparrow}$ can be felt far downwind for power-law like $P_0(\ell)$.

Panel (c) of Figure S2 shows how Gaussian/delta-like $P_0$ effects ripple growth. The growth rates for the salton-repton model for the Gaussian-like $P_0$ gives rise to a fastest growing wavelength that scales with the average hoplength $\lambda \sim \langle\ell\rangle$ and approaches $\lambda \to 4\langle\ell\rangle$ as $P_0 \to \delta(\ell - \langle\ell\rangle)$ (Eq. S24). However for the more realistic power-law like $P_0$ there is no maxima and thus no wavelength selection (the smallest wavelength grows the fastest). Also highlighted in Figure S2c (insets) is the fact that in the local regime $(\ell/\lambda \ll 1)$ the form of $P_0$ has little importance, however in the nonlocal regime $(\ell/\lambda \gtrsim 1)$ the form of $P_0$ is crucial: For $(\ell/\lambda \ll 1)$ grains are redistributed locally on the ripple face—thus the flux $Q$ is lagged from $\varphi_{\uparrow}$ by a small amount (set by $(\ell/\lambda)$) relative to $\lambda$. However, in the $(\ell/\lambda \gtrsim 1)$ regime, $Q(x)$ at any spot $x$ is nonlocally related to $\varphi_{\uparrow}(x)$ far upwind. Thus for the unrealistic $P_0(\ell)$, the peak of the Gaussian-like density of grains tends to cause a large lag between $\varphi_{\uparrow}$ and $Q$, shifting the maximum of $Q$ past the crest and leading to erosion of the crest and stabilization of small wavelengths. This can also be viewed as grains filling up the downwind troughs. However for the realistic $P_0(\ell)$, the flux relaxes so slowly to changes in $\varphi_{\uparrow}$ (Fig. S2b), because the grains are weighted heavily at $\ell = 0$, that the flux maximum never passes the crest, even as $\lambda \to 0$. Therefore, the salton-repton model, without inclusion of an impact-ejection lag, does not include a stabilizing mechanism and
therefore cannot explain wavelength selection for impact ripples.

## E. Modulation Rate Approximations

To gain insight into the behavior of the flux modulations $\mathcal{F}_\uparrow$ and $\mathcal{A}(\ell)$ we look at the dynamics around the marginal stability point $\Omega(k_*) = 0$ where the bed perturbation $k_*$ does not grow or propagate. This occurs when the magnitude and phase of the vertical fluxes match, $\mathcal{F}_\downarrow = \mathcal{F}_\uparrow$, which occurs for wavelengths of order $\langle\ell\rangle$ where the hoplength amplification effects for $\mathcal{A}(\ell \sim \lambda)$ are significant. Note that $k_*$ is of order (around 2 times) the fastest growing mode and thus the dynamics around $k_*$ tell us about the dynamics that give rise to the ripples that emerge from a flat bed.

Using Equation (S18) we conclude that the marginal vertical flux modulations are $\mathcal{F}_* = i\langle\mathcal{L}\rangle/(\tilde{P}_0 - 1)$. If we approximate the hoplength modulation as being purely geometric (Eq. (S17)) we can approximate the flux modulation around the marginal stability mode as

$$ \mathcal{F}_\uparrow \sim \mathcal{F}_* + \delta\mathcal{F}_\uparrow \sim i\langle\cot\theta\rangle + \delta\mathcal{F}_\uparrow $$

where for illustration we take $\mathcal{F}_* = i\langle\cot\theta\rangle$ and absorb the error in $\delta\mathcal{F}_\uparrow$. We can see that the correction to the vertical flux $\delta\mathcal{F}_\uparrow$ is responsible for ripple growth around $k_*$:

$$ \frac{\Omega}{k\varphi_0} = (\tilde{P}_0 - 1)(\mathcal{F}_\uparrow - \mathcal{F}_*) = (\tilde{P}_0 - 1)\delta\mathcal{F}_\uparrow. $$

The point here is that around the mode $k_*$ (of order the fastest growing mode), the hoplength amplification forcing is strong and thus tends to maximize the vertical flux $\varphi_\uparrow$ at the upwind inflection point (to 0-th order). We can see the amplified forcing around $k_*$ by observing $\mathcal{A}_*(\ell|\mathcal{F}_*)$ using equation (S21) for $\ell/\lambda \gtrsim 1$:

$$ \frac{\mathcal{A}_*}{\langle\cot\theta\rangle} \sim i(e^{-ik\ell} + 1) $$

where we have ignored the small contributions from $-k^{-1}P'_0/P_0 \sim (k\ell)^{-1}$. This expression highlights that $\mathcal{A} = \mathcal{A}_* + \delta\mathcal{A}$ is approximately a unique function of $\ell/\lambda$, thus showing that the forcing from trajectories $\ell \sim \lambda$ is approximately scale invariant, $\mathcal{A}(\ell \sim \lambda) \sim \mathcal{A}_*(\ell/\lambda)$. The forcing amplification is seen via $\mathcal{A}_*(\ell = \lambda) \simeq 2i\langle\cot\theta\rangle$ which is maximized around the upwind inflection point ($\pi/2$ out of phase; $i = e^{i\pi/2}$). This scale-invariant forcing is observed in DEM measurements of $\mathcal{A}(\ell)$ seen in Methods Figure M3.

## F. Multiple-scales: Impact-Ejection Length Selects the Wavelength

Because we are only interested in the regime where $k\epsilon \ll 1$ (namely, $\lambda \gg d$) we can approximate $\tilde{n}_\epsilon(\ell) = \int n_\epsilon(\epsilon|\ell)e^{-ik\epsilon}d\epsilon \sim n_\uparrow(\ell)(1 - ik\epsilon(\ell))$. The modulation rate $\mathcal{F}_\uparrow$ can then be expanded into a large-scale/0-th or-
<figure><img src="image_11.png" /><figcaption>FIG. S3. Magnitude and phase of the ejection flux modulation rate $\mathcal{F}_{\uparrow} = \hat{\phi}_{\uparrow}/(\phi_0 k \hat{Z})$. Symbols are measurements from DEM simulations and lines are approximations using geometrical hoplength changes $\mathcal{L} = \cot\theta(l)(e^{-ikl} - 1)$ in Equation S22. The dashed lines show the large-scale/0-th order modulation $\mathcal{F}_{\uparrow}(\epsilon = 0) \simeq f_{\uparrow}$ (Eq. S25). The disagreement between measured and approximated $\mathcal{F}_{\uparrow}$ is likely due to measurement error and missing contributions to $\mathcal{L}$ (see Methods).</figcaption></figure>

der modulation, which reflects the $\Delta\ell$ modulations of the impact forcing, and a small-scale/1-st order modulation contribution reflecting the downwind shift in ejections from $\epsilon$. We can approximate the modulation rate generally as

$$ \mathcal{F}_{\uparrow} \approx f_{\uparrow}(1 - ik\langle\epsilon_{\Delta\ell}\rangle_n) \approx f_{\uparrow}e^{-ik\langle\epsilon_{\Delta\ell}\rangle_n} \quad (S25) $$

where $f_{\uparrow} = \rho_n^{-1}(-i\langle\mathcal{L}\rangle_n + k^{-1}\langle n'_{\uparrow}\mathcal{L}\rangle)$ is the large-scale/0-th order modulation which is an approximate reflection of the forcing $f_{\uparrow} \sim \int n_{\uparrow}(l)A(l|f_{\uparrow})P_0(l)d\ell$—or put another way, $f_{\uparrow} \sim \mathcal{F}_{\uparrow}(\epsilon = 0)$. Taking into account the nonlocal impact-ejection lag $\epsilon$, $f_{\uparrow}$ then becomes phase-shifted downwind by the small length scale $\langle\epsilon_{\Delta\ell}\rangle_n = \langle\bar{\epsilon}\mathcal{L}\rangle_n/\langle\mathcal{L}\rangle_n$ which is the average of the lag $\bar{\epsilon}$ weighted by the changes in hoplengths over the distribution $P_n = n_{\uparrow}P_0$ where we define $\langle\cdot\rangle_n = \int(\cdot)P_n d\ell$.
1. Local Limit: $\lambda \gg \langle \ell \rangle$

In the local regime ($\lambda \gg \langle \ell \rangle$) the geometrical hoplength modulation becomes determined by the local slope such that $\mathcal{L} \sim -ik\ell \cot \theta$ and effects from $\epsilon$ are minor. Therefore the vertical flux modulation in this limit behaves as (Fig. S3)

$$ \lim_{k\langle\ell\rangle\to 0} \frac{\mathcal{F}_{\uparrow}}{\langle\cot\theta_{\ell}\rangle_n} \sim i - k^{-1}\langle n'_{\uparrow\ell\theta} \rangle $$

where $\langle \cot \theta_{\ell} \rangle_n = \langle \ell \cot \theta \rangle_n / \langle \ell \rangle_n$ and $\langle n'_{\uparrow\ell\theta} \rangle = \langle n'_{\uparrow\ell} \ell \cot \theta \rangle / \langle \ell \cot \theta \rangle_n$. This shows that the vertical flux tends to increase with the slope $\mathcal{F}_{\uparrow} \sim i$ as in the idealized salton-repton model, but for low slopes (large $\lambda$) the vertical flux modulation also tends to be maximized in the ripples trough do to the gradients in ejections $n'_{\uparrow}(\ell)$—that become important when the impacts are shorted or lengthened (impact energy lost or gained) from local slope effects.

Because the effects of $\epsilon$ are vanishing in the limit as $\langle \ell \rangle / \lambda \to 0$, ripples grow and propagate according to

$$ \begin{gathered} \lim_{k\langle\ell\rangle\to 0} \frac{\sigma}{k^2 Q_0} \sim \langle\cot\theta_{\ell}\rangle_n - \langle\cot\theta_{\ell}\rangle > 0 \\ \lim_{k\langle\ell\rangle\to 0} \frac{c}{\varphi_0} \sim -\frac{\langle\ell\rangle}{\langle\ell\rangle_n} \langle n'_{\uparrow}\ell \cot\theta\rangle < 0, \end{gathered} \qquad (S26) $$

showing that flat bed ripple growth for large wavelengths occurs through antidiffusion ($\partial_t Z \sim -Q_0 \partial_{xx} Z$) at a rate given by the flux $Q_0 = \langle \ell \rangle \varphi_0$ and that the speeds in this regime are purely negative (antiripples).

2. Nonlocal Limit: $\lambda \ll \langle \ell \rangle$

For wavelengths around the fastest growing wavelength, the flux response is nonlocal ($\lambda \lesssim \langle \ell \rangle$; Fig. 4, S2) and thus the impact flux amplification/focusing on the upwind inflection point dominates the forcing. This renders the large scale modulation of the ejection flux (from $\Delta\ell$ effects) maximized around the inflection point, $f_{\uparrow} \sim i \langle \cot \theta \rangle_n$ (Fig. S3). Thus, in the far nonlocal regime around the maximum of $\sigma$ we have

$$ \begin{gathered} \lim_{k\langle\ell\rangle\to\infty} \frac{\sigma}{k\varphi_0\langle\cot\theta\rangle_n} \sim \rho_{\text{Im}} + \frac{\langle\mathcal{L}\rangle_{\text{Im}}}{\langle\cot\theta\rangle_n} - \rho_{\text{Re}}k\langle\varepsilon_{\theta}\rangle_n \\ \lim_{k\langle\ell\rangle\to\infty} \frac{c}{\varphi_0\langle\cot\theta\rangle_n} \sim \rho_{\text{Re}} + \frac{\langle\mathcal{L}\rangle_{\text{Re}}}{\langle\cot\theta\rangle_n} + \rho_{\text{Im}}k\langle\varepsilon_{\theta}\rangle_n \end{gathered} \qquad (S27) $$

where $\rho_{\text{Re}} = 1 - \tilde{P}_0|_{\text{Re}} < 1$ and $\rho_{\text{Im}} = -\tilde{P}_0|_{\text{Im}} \ll 1$, are both positive and weak functions of $k\langle\ell\rangle$ (for $\lambda \ll \langle \ell \rangle$), and $\langle \varepsilon_{\theta} \rangle_n = \langle \bar{\varepsilon} \cot \theta \rangle_n / \langle \cot \theta \rangle_n$. This recovers the empirical scaling found by Duran et al. (2014) [9] with $\sigma \propto Ak - Bk^2 \langle \varepsilon_{\theta} \rangle_n$. Here $A = \langle \cot \theta \rangle_n \rho_{\text{Im}} + \langle \mathcal{L} \rangle_{\text{Im}}$ and $B = \langle \cot \theta \rangle_n \rho_{\text{Re}}$ are weak functions of $k$ in the nonlocal regime.
This shows that we can approximate that the fastest growing wavelength scales roughly with the average impact-ejection lag length, $\lambda \sim \frac{B}{A}\langle\epsilon_{\theta}\rangle_n$. Inserting this scaling into $c(k)$ we find the speed scales with transport dynamics through $\varphi_0$ and $\langle\ell\rangle$, $c \sim \varphi_0 f(\langle\ell\rangle)$. Thus in this regime the wavelengths are selected by the granular dynamics scale that depends very little on transport while the associated propagation speeds depend very strongly on the transport regime. Note that these approximations are most accurate for exceedingly large $\langle\ell\rangle$ however Equations (S27) are still a reasonable approximation for Mars and Earth-like regimes (Figure 4) and support the general scaling of the model (Methods Figure M4). Interestingly, the term $\rho_{Re} + \langle\mathcal{L}\rangle_{Re}/\langle\cot\theta\rangle_n$ in $c(k)$ in eq. (S27) is positive for $\lambda \ll \langle\ell\rangle$ and negative for $\lambda \gtrsim \langle\ell\rangle$. Thus the approximation for $c$ in (S27) can approximate the ripple-antiripple scaling transition.

### G. Dispersion Relation in Terms of The Scale-Invariant Amplification Rate

Another more elegant way of seeing the wavelength selection mechanisms is by using the growth rate formalism in terms of $\mathcal{A}$ (Eq. S20). Like in the Methods, we define $\mathcal{A} - \delta\mathcal{A} = \mathcal{R} + i\mathcal{S}$ where $\mathcal{R}$ and $\mathcal{S}$ are the in-phase and out-of-phase contributions to the impact flux amplification, respectively, determined by hoplength changes $\mathcal{L}(k\ell)$. Here $\delta\mathcal{A} \simeq -ik\langle\epsilon_{\Delta\ell}\rangle_n f_{\uparrow}e^{-ik\ell}$ is the small perturbation to the nearly scale-free forms of $\mathcal{R}(k\ell)$ and $\mathcal{S}(k\ell)$ (Fig. M3). We note that in the Methods we ignore the small effects of $\delta\mathcal{A}$ for brevity and subtract $\delta\mathcal{A} \simeq -ik\langle\epsilon_{\Delta\ell}\rangle_n f_{\uparrow}e^{-ik\ell}$ from direct DEM measurements of $\mathcal{A}$ in Methods Figure M3 to obtain the DEM values of $\mathcal{R}$ and $\mathcal{S}$ presented. Using these definitions we can express the dispersion relation as

$$ \begin{aligned} \frac{\sigma}{k\varphi_0} &= \langle(1-n_{\uparrow})\mathcal{R}\rangle - k\langle\bar{\epsilon}n_{\uparrow}\mathcal{S}\rangle \\ \frac{c}{\varphi_0} &= \langle(n_{\uparrow}-1)\mathcal{S}\rangle - k\langle\bar{\epsilon}n_{\uparrow}\mathcal{R}\rangle \end{aligned} \qquad (S28) $$

ignoring the contributions from $\delta\mathcal{A}$ for conceptual clarity. We saw in SI Section II E for large $\ell/\lambda \gtrsim 1$ around $k \sim k_*$ (small $\lambda/\langle\ell\rangle$), $\mathcal{A} \sim \mathcal{A}_*$ such that $\mathcal{R}_* \propto \sin k\ell$ and $\mathcal{S}_* \propto 1 + \cos k\ell$ (Fig. M3). Thus for nonlocal $\lambda/\langle\ell\rangle$, because most of the grains are leaving from the upwind face ($f_{\uparrow} \propto i$), $\mathcal{R}$ is peaked for small hoplengths around $\ell \sim \lambda/4$, acting like a 'repton' amplification. Similarly, in the nonlocal regime the focusing effects ensure that $\mathcal{S}$ is peaked around $\ell \sim \lambda$, acting like a 'salton' amplification.

This way of stating the dispersion relation ultimately results in the same expressions derived above but now has a particularly intuitive meaning: The destabilizing mechanism now represents the average number of grains that remain on the crest after traveling from the upwind face $\langle(1-n_{\uparrow})\mathcal{R}\rangle$. This effect is countered by the stabilization mechanism from the increase number of grains being eroded from the crest due to the down wind shift in the ejection flux from $k\langle\bar{\epsilon}n_{\uparrow}\mathcal{S}\rangle \sim k\langle\epsilon_{\Delta\ell}\rangle_n\langle n_{\uparrow}\mathcal{S}\rangle$—this effect increases with the number of grains leaving the upwind face, $\langle n_{\uparrow}\mathcal{S}\rangle$.

Because the number of ejections/rebounds $n_{\uparrow}(\ell)$ is small for small hoplengths (Fig. 3), the speed contribution from $k\langle\bar{\epsilon}n_{\uparrow}\mathcal{R}\rangle$ is negligible thus leaving the propagation speed to be determined by the number of grains eroded from or deposited on the upwind face $\langle(n_{\uparrow}-1)\mathcal{S}\rangle$ which is positive in the nonlocal regime where $f_{\uparrow} \propto i$ because grains are
being distributed from the upwind face over to the back of the ripple. However, if the wavelength is selected in the local regime (small $\langle l \rangle$) then the speeds can be negative.

With these considerations we arrive at the approximation

$$ \frac{\sigma}{k\varphi_0} \sim \langle (1 - n_\uparrow)\mathcal{R} \rangle - k\langle \varepsilon_{\Delta l} \rangle_n \langle n_\uparrow\mathcal{S} \rangle \\ \frac{c}{\varphi_0} \sim \langle (n_\uparrow - 1)\mathcal{S} \rangle \qquad (S29) $$

which reduces to the main text finding that

$$ \frac{\lambda}{\langle \varepsilon_{\Delta l} \rangle_n} \propto \frac{\langle n_\uparrow \mathcal{S} \rangle}{\langle (1 - n_\uparrow) \mathcal{R} \rangle} $$

with a proportionality factor $a(\langle l \rangle)$ associated with the contribution from $\delta A$.

## H. The Transport Regime Sets the Propagation Speed

As we have seen above the fastest growing wavelength can be written generally as $\lambda/\langle \varepsilon_{\Delta l} \rangle_n \simeq f(\langle l \rangle/d)$. The factor $f(\langle l \rangle/d)$ is an almost unique function of $\langle l \rangle/d$ and so is $\langle \varepsilon_{\Delta l} \rangle_n/d$. Thus the dispersion relation can be expressed as a function of two groupings, $k\langle l \rangle$ and $\langle l \rangle/d$:

$$ \frac{\Omega}{k\varphi_0} = G(k\langle l \rangle, \langle l \rangle/d). $$

The fastest growing wavelength follows from the condition $\partial_k \sigma = 0$ which renders $k\langle l \rangle = g(\langle l \rangle/d)$ and thus the speed is found to be a unique function of $\langle l \rangle/d$:

$$ \frac{c}{\varphi_0} = F(\langle l \rangle/d) $$

recovering the previous approximation that while $\lambda$ scales with a granular length $\langle \varepsilon_{\Delta l} \rangle_n \propto d$, the speed $c$ scales with the transport regime, $\varphi_0$ and $\langle l \rangle/d$ (Fig. M4).