Unified threshold model suggests sediment transport by wind on Triton, Pluto, and comet 67P
Thomas P\"ahtz
$^{1
2
When Voyager~2 imaged the surface of Neptune's moon Triton in 1989, it revealed the occurrence of surface streaks that are possibly of aeolian origin (i.e., wind-formed) \cite{Hansenetal90,SaganChyba90. Likewise, New Horizons imaged surface features that have been tentatively interpreted as possible wind streaks when it passed Pluto in 2015 \cite{Sternetal15. Moreover, Rosetta imaged what looked like aeolian ripples and dunes on the comet 67P/Churyumov-Gerasimenko (67P) in 2014 \cite{Thomasetal15a,Thomasetal15b. However, whether these surface features formed due to aeolian sediment transport remains a mystery \cite{SaganChyba90,Sternetal15,Thomasetal15a,Thomasetal15b because the atmospheres on these planetary bodies are extremely thin. In fact, it has been estimated that average 1m winds of more than 500km/h are required to lift sediment directly from the surface on Triton and Pluto \cite{Mooreetal15, where winds are weaker than on Earth \cite{SaganChyba90,Gladstoneetal16. Here, using physical modelling, we drastically lower these estimates. We predict that sediment transport can be sustained under winds that are weaker than the strongest estimated winds occurring on these planetary bodies. The main reason is entrainment of sediment from the surface through impacts of transported particles, which has already been described as the reason for low thresholds on Mars \cite{Kok10a. This mechanism requires that a few particles initially enter the air flow through a process different from wind erosion. Geyser-like eruptions on Triton and airfall of particles orbiting comet 67P, reported by the Voyager~2 and Rosetta flyby, respectively, provide the air flow with plenty of particles to initiate transport. For Pluto, impacts of micrometeoroids is the most probable initiation mechanism. We estimate the recurrence time of relevant impacts is several orders of magnitude smaller than the timescale associated with the formation of large-scale landforms in the absence of aeolian and fluvial processes. We thus conclude surface features on Triton, Pluto, and 67P very likely formed due to aeolian sediment transport.
Sedimentthreshold.tex
%% Template for a preprint Letter or Article for submission
%% to the journal Nature.
%% Written by Peter Czoschke, 26 February 2004
%%

\documentclass{nature}

\pdfoutput=1
\usepackage{graphicx}
\usepackage{epstopdf}
\usepackage{amssymb}
\usepackage{mathtools}
\usepackage{xfrac}
\usepackage{color}

\def \d{\mathrm{d}}

\definecolor{red}{rgb}{1,0,0} 
\newcommand{\ed}[1]{ \textcolor{red}{#1} }
\newcommand{\cred}{\color{red}}
\newcommand{\cblack}{\color{black}}

%% make sure you have the nature.cls and naturemag.bst files where
%% LaTeX can find them

%\bibliographystyle{naturemag}

\title{Unified threshold model suggests sediment transport by wind on Triton, Pluto, and comet 67P}

%% Notice placement of commas and superscripts and use of &
%% in the author list

\author{Thomas P\"ahtz,$^{1,2}$ Orencio Dur\'an$^3$}


\begin{document}

\maketitle

\begin{affiliations}
 \item Institute of Port, Coastal and Offshore Engineering, Ocean College, Zhejiang University, 866 Yu Hang Tang Road, 310058 Hangzhou, China
 \item State Key Laboratory of Satellite Ocean Environment Dynamics, Second Institute of Oceanography, 36 North Baochu Road, 310012 Hangzhou, China
 \item Department of Physical Sciences, Virginia Institute of Marine Sciences, College of William and Mary, PO Box 1386, Gloucester Point, 23062 Williamsburg, USA
\end{affiliations}

\begin{abstract}
When Voyager~2 imaged the surface of Neptune's moon Triton in 1989, it revealed the occurrence of surface streaks that are possibly of aeolian origin (i.e., wind-formed) \cite{Hansenetal90,SaganChyba90}. Likewise, New Horizons imaged surface features that have been tentatively interpreted as possible wind streaks when it passed Pluto in 2015 \cite{Sternetal15}. Moreover, Rosetta imaged what looked like aeolian ripples and dunes on the comet 67P/Churyumov-Gerasimenko (67P) in 2014 \cite{Thomasetal15a,Thomasetal15b}. However, whether these surface features formed due to aeolian sediment transport remains a mystery \cite{SaganChyba90,Sternetal15,Thomasetal15a,Thomasetal15b} because the atmospheres on these planetary bodies are extremely thin. In fact, it has been estimated that average 1m winds of more than 500km/h are required to lift sediment directly from the surface on Triton and Pluto \cite{Mooreetal15}, where winds are weaker than on Earth \cite{SaganChyba90,Gladstoneetal16}. Here, using physical modelling, we drastically lower these estimates. We predict that sediment transport can be sustained under winds that are weaker than the strongest estimated winds occurring on these planetary bodies. The main reason is entrainment of sediment from the surface through impacts of transported particles, which has already been described as the reason for low thresholds on Mars \cite{Kok10a}. This mechanism requires that a few particles initially enter the air flow through a process different from wind erosion. Geyser-like eruptions on Triton and airfall of particles orbiting comet 67P, reported by the Voyager~2 and Rosetta flyby, respectively, provide the air flow with plenty of particles to initiate transport. For Pluto, impacts of micrometeoroids is the most probable initiation mechanism. We estimate the recurrence time of relevant impacts is several orders of magnitude smaller than the timescale associated with the formation of large-scale landforms in the absence of aeolian and fluvial processes. We thus conclude surface features on Triton, Pluto, and 67P very likely formed due to aeolian sediment transport.
\end{abstract}

Although previous studies have proposed analytical models that predict planetary transport cessation thresholds \cite{ClaudinAndreotti06,Kok10b,Paehtzetal12}, their predictions can differ by up to one order of magnitude when applied to different planetary conditions (Extended Data Figure~1). The main reason is that they incorporate different assumptions of how the entrainment of sediment through impacts of transported particles affect the average particle dynamics. None of these assumptions has been fully justified because experiments today are far too imprecise to provide the needed information. Using numerical simulations of sediment transport with a Discrete Element Method model coupled with the Reynolds Averaged Navier-Stokes equations \cite{Duranetal12}, we have recently overcome these difficulties and provided a simple physical parametrization of the effects of particle-bed impcats \cite{PaehtzDuran16a}. Based on this result, we here propose a unified analytical model of the cessation threshold of sediment transport that can be applied to arbitrary environments and thus provides estimates of the minimal fluid speeds required to shape fluid-sheared surfaces as diverse as sea floors, riverscapes, and wind-blown loose planetary surfaces (Fig.~\ref{Fluid_sheared_surfaces}), including those on Triton, Pluto, and comet 67P. A major strength of the model, in contrast to previous ones, is that it reproduces measurements in both Earth's atmosphere and under water (Fig.~\ref{Modelvsmeasurements}). This is, indeed, a requirement for any model based on the impact entrainment mechanism since recently reported evidence by three independent studies strongly indicates that turbulent subaqueous sediment transport is also sustained through particle-bed impacts \cite{PaehtzDuran16a,Clarketal15,Maurinetal15}, contradicting the previous state of the art. This evidence consists of direct visualization of impact entrainment events \cite{PaehtzDuran16a} (see also supplementary movies), a constant average particle slip velocity \cite{PaehtzDuran16a}, and a transport cessation threshold that is significantly below the initiation threshold \cite{Clarketal15,Maurinetal15}. The latter two are widely recognized as indirect evidence of predominant impact entrainment \cite{Koketal12}.

%
\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Fig1.png}
 \end{center}
 \caption{{\bf Example of bedforms for different planetary conditions}. (a) under water, (b) Earth, (c) Mars, (d) possibly on Pluto, and (e) possibly on comet 67P/Churyumov-Gerasimenko. Images are from (a) https://en.wikipedia.org/wiki/Giant\_current\_ripples (license CC BY-SA 3.0), (b) https://en.wikipedia.org/wiki/Desert (license CC BY-SA 3.0), (c) NASA/JPL-Caltech/University of Arizona (Id: ESP\_011909), (d) NASA/JHUAPL/SWRI (Id: Snakeskin), and (e) ESA/Rosetta/MPS (Id: 333962).}
 \label{Fluid_sheared_surfaces}
\end{figure}
%

In general, sediment transport in such a variety of conditions depends on the fluid shear velocity $u_\ast$ or shear stress $\tau = \rho_f u^2_\ast$, gravitational constant $g$, particle density $\rho_p$, fluid density $\rho_f$, kinematic fluid viscosity $\nu$, mean particle diameter $d$, mean free path of the fluid molecules $\lambda$, and surface tension $\gamma$ and Young modulus $E$ of the particles. From dimensional analysis, the rescaled shear stress, or Shields parameter $\Theta=u_\ast^2/[(s-1)gd]$, at the transport threshold $\Theta^{\mathrm{ex}}_t=u^{\mathrm{ex}2}_t/[(s-1)gd]$ (obtained from extrapolation to vanishing transport, see Methods) is then a function of the particle-fluid-density ratio $s=\rho_p/\rho_f$, particle Reynolds number $\mathrm{Re}_p=\sqrt{(s-1)gd^3}/\nu$, cohesion number $C=d^{-1}\gamma^{3/5}E^{-1/5}[(\rho_p-\rho_f)g]^{-2/5}$, and Knudsen number $\mathrm{Kn}=\lambda/d$.

In contrast to the onset of sediment transport, where particles at rest are set in motion only by fluid forces, when approaching transport cessation, particles at the bed enter the flow mainly due to the impact of moving particles, which requires lower flow velocities. Equilibrium is reached when, on average, only one particle per impact enters the flow. The transport cessation threshold is then defined as the limit of vanishing number of transported particles per unit area \cite{Duranetal12} (Methods). However, in this extrapolated state, particles are still in motion, meaning that physical quantities characterizing particle trajectories do not vanish, and there is a well defined equilibrium transport layer. The threshold function $\Theta^{\mathrm{ex}}_t(s,\mathrm{Re}_p,C,\mathrm{Kn})$ is obtained from a set of equations describing the physical properties of the equilibrium transport layer close to the cessation threshold. These equations relate four dimensionless variables characterizing the transport layer: the average horizontal and vertical particle velocity ($V_x=\overline{v_x}/\sqrt{(s-1)gd}$ and $V_z=\sqrt{\overline{v_z^2}}/\sqrt{(s-1)gd}$, respectively), the transport layer height $Z=(\overline{z}-z_s)/d$ relative to the bed surface $z_s$, and the average fluid velocity $U_x=\overline{u_x}/\sqrt{(s-1)gd}$. The overbar denotes an average over all particles in the transport layer (i.e., located above $z_s$). A summary of all dimensionless quantities and their meaning is shown in Table~\ref{DimensionlessQuantities}. A general description of the model is presented below, see Methods for the details.

\begin{table}
\resizebox*{\textwidth}{!}{%
    \begin{tabular}{| p{0.4\columnwidth} | p{0.6\columnwidth} |}
    \hline
     Dimensionless quantity & Physical meaning \\ \hline
		 Threshold Shields parameter $\Theta^{\mathrm{ex}}_t$ & Measure of the importance of drag forces at the cessation threshold relative to the submerged weight of particles. \\ \hline
		 Particle-fluid-density ratio $s$ & Main parameter distinguishing planetary environments. \\ \hline
		 Particle Reynolds number $\mathrm{Re}_p$ & Inverse measure of the importance of viscous drag forces relative to the submerged weight of particles. Small values are typical for rarefied environments. \\ \hline
		 Dimensionless viscous length $Z_\nu \equiv 1/(\mathrm{Re}_p\sqrt{\Theta^{\mathrm{ex}}_t})$ & Measure of the size of the viscous sublayer of the turbulent boundary layer at the cessation threshold relative to particle size. \\ \hline
		 Cohesion number $C$ & Measure of the importance of cohesive interparticle forces. Only important for small particles. \\ \hline
		 Knudsen number $\mathrm{Kn}$ & Mean free path of fluid molecules relative to particle size. Only important in rarefied environments. \\ \hline
		 Dimensionless average horizontal and vertical particle velocity $(V_x,V_z)$ & Characterizes particle motion near threshold conditions. \\ \hline
		 Dimensionless average horizontal fluid velocity $U_x$ & Measure of the typical fluid velocity in the transport layer near threshold conditions. \\ \hline
		 Dimensionless transport layer height $Z$ & Measure of the typical elevation of particle trajectories. \\	 
    \hline
    \end{tabular}
}
		\caption{Physical meaning of dimensionless parameters and model quantities.}
		\label{DimensionlessQuantities}
\end{table}

Physical arguments combined with direct simulations of cohesionless sediment transport \cite{Duranetal12} suggest the equilibrium state close to the cessation threshold has several characteristics. First, the `friction coefficient' (the ratio of granular shear stress and pressure, $\mu$) is relatively constant at the interface between the dynamic transport layer and the quasi-static sediment bed \cite{Bagnold56,Duranetal12,Paehtzetal12} ($\mu_b=\mu(z_s)=\mathrm{const}$) (Extended Data Figure~2). The friction coefficient $\mu_b$ describes the particle resistance to shear due to intergranular forces and collisions with the bed. Because this resistance increases with cohesive forces, $\mu_b$ depends on the cohesion number $C$ (Methods). From the momentum balance within the transport layer, a constant friction coefficient means that the average horizontal fluid drag acceleration $\overline{a^d_x}$ is proportional to the vertical buoyancy-reduced gravitational acceleration $\tilde g=(1-1/s)g$ \cite{Paehtzetal12}. For single particles in non-rarefied systems, $a^d_x/\tilde g$ is a well-known function of $\mathrm{Re}_p$ and $U_x-V_x$ \cite{Camenen07}. However, averaging and rarefaction effects introduce additional dependencies on $V_z$ and $\mathrm{Kn}$ \cite{Croweetal12}, respectively. Thus $\mu$ has the functional ($f_1$) form (Methods), 
%
\begin{equation}
 \mu_b(C) = \overline{a^d_x}/\tilde g = f_1(\mathrm{Re}_p,U_x-V_x,V_z,\mathrm{Kn}). \label{mub1} 
\end{equation}
%
Notice that it is crucial to include rarefaction (i.e., large $\mathrm{Kn}$) effects on the fluid drag as they can modify $\overline{a^d_x}$ by up to an order of magnitude on Triton and Pluto and several orders of magnitude on comet 67P for typical conditions. Due to their thin atmospheres, the air molecules have a mean free path that is comparable (Triton and Pluto) or much larger (67P) than the particle size. 

Second, close to the cessation threshold, the rescaled fluid velocity above the bed ($u_x(z)/u^{\mathrm{ex}}_t$) is mostly undisturbed by particle motion and approaches a well-known function (``the law of the wall'') of the dimensionless elevation $z/d$ rescaled by the viscous length $Z_\nu = 1/(\mathrm{Re}_p\sqrt{\Theta^{\mathrm{ex}}_t})$. The average rescaled fluid velocity within the transport layer ($\overline{u_x}/u^{\mathrm{ex}}_t=U_x/\sqrt{\Theta^{\mathrm{ex}}_t}$) can then be approximated by the fluid velocity evaluated at the transport layer height $Z$
%
\begin{equation}
 	U_x = \sqrt{\Theta^{\mathrm{ex}}_t}f_2(Z/Z_\nu), \label{MeanU}
\end{equation}
where the ``the law of the wall'' velocity profile $f_2$ is defined in Methods. 

Third, the scaling of the average particle velocity $V_x$ depends on the size of the transport layer. When the transport layer is within the viscous sublayer of the turbulent boundary layer ($Z/Z_\nu \lessapprox 7$), $V_x$ scales with the characteristic fluid velocity within the viscous sublayer: $V_x \propto U_x \approx \sqrt{\Theta^{\mathrm{ex}}_t} Z/Z_\nu$ (Extended Data Figure~3). When the transport layer extends beyond the viscous sublayer, the scale of the particle velocity is dominated by the characteristic fluid velocity in the logarithmic region of the velocity profile: $V_x \propto \sqrt{\Theta^{\mathrm{ex}}_t}$. The transition between these two regimes is captured by a phenomenological expression obtained from transport simulations (Extended Data Figure~3b),
%
\begin{equation}
	V_x/\sqrt{\Theta^{\mathrm{ex}}_t}=\beta_1\sqrt{1-\exp\left[-\beta_2\left(U_x/\sqrt{\Theta^{\mathrm{ex}}_t}\right)^2\right]}, \label{VU}
\end{equation}
%
where $\beta_1=4.8$ and $\beta_2=0.027$ are fitting parameters.

\cblack

Fourth, from particles' vertical momentum balance, the transport layer height $Z$ increases with the granular temperature $V_z^2$, with a small offset $Z_c$ due to interparticle contacts near the surface (Methods),
%
\begin{equation}
	Z = Z_c + s V_z^2. \label{Meanz}
\end{equation}
%

Finally, from the definition of the granular stress tensor and the constancy of the friction coefficient at the bed, the vertical and horizontal average particle velocities are found to be proportional to each other (Methods), 
%
\begin{equation}
	V_z = \alpha V_x /\mu_b(C). \label{TzV} 
\end{equation}
%

Equations (\ref{Meanz}) and (\ref{TzV}) are consistent with data from cohesionless, non-rarefied ($C=\mathrm{Kn}=0$) transport simulations (Extended Data Figures~4 and 5), from which $Z_c = 0.9\mu_b^{-1}\Theta$ and $\alpha = 0.18$ are obtained.

%
\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=0.7\columnwidth]{Fig2.pdf}
 \end{center}
 \caption{{\bf Validation of model predictions.} Threshold Shields parameter $\Theta^{\mathrm{ex}}_t$ versus the particle Reynolds number $\mathrm{Re}_p$ predicted by the analytical model (solid lines) and measured in water \cite{Buffingtonmontgomery97,Ouremietal07} (blue, open symbols, dashed line) and Earth's atmosphere \cite{Bagnold37,Chepil45,Creysselsetal09,Hoetal11} (brown, open symbols). Simulations by the cohesionless transport model~\cite{Duranetal12} are also shown (filled symbols). The review of subaqueous threshold data by Ref.~\cite{Buffingtonmontgomery97} reports cessation threshold values (measured via extrapolation) of numerous studies, we thus choose to represent the average and standard error ($95\%$ confidence) over all available measurements (blue symbols). For water/water-oil mixtures, cohesionless predictions were used because measurements correspond to very-weak-cohesion conditions ($d\geq132\mu$m) as cohesive forces under water are much weaker than in dry environments due to the lack of electrostatic forces \cite{Koketal12}. The model overestimates the constant value of the viscous subaqueous threshold (measured via extrapolation \cite{Ouremietal07}), which is expected because hindrance effects are neglected (supplementary online text). For air, Refs.~\cite{Bagnold37} and \cite{Chepil45} directly measured the threshold as the smallest values of $\Theta$ below which sediment transport ceases, whereas for Refs.~\cite{Creysselsetal09,Hoetal11}, we extrapolated the threshold from measurements of the sediment transport rate.}
 \label{Modelvsmeasurements}
\end{figure}
%

%
\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Fig3.pdf}
 \end{center}
 \caption{{\bf Transport threshold predictions for planetary bodies.} Predicted threshold shear velocity $u^\mathrm{ex}_t$ as a function of the mean particle diameter $d$ for various planetary environments (See Methods for details.) Symbols correspond to measurements in Earth's atmosphere \cite{Bagnold37,Chepil45,Creysselsetal09,Hoetal11} and the shaded areas to regions of possible sediment transport inferred from measurements/estimations of the strongest winds occurring on the respective planetary body (Methods).}
 \label{Planetpredictions}
\end{figure}
%

The analytical model (Eqs.~(\ref{mub1}-\ref{TzV})) predicts the whole range of cessation thresholds obtained from cohesionless, non-rarefied transport simulations (Extended Data Figure~6). It also reproduces measurements in such different conditions as water/water-oil mixtures and Earth's atmosphere \cite{Bagnold37,Chepil45,Buffingtonmontgomery97,Ouremietal07,Creysselsetal09,Hoetal11} without fitting any model parameter to the threshold data (Fig.~\ref{Modelvsmeasurements}).

When applied to various planetary  environments (Fig.~\ref{Planetpredictions}), we find the predictions for Triton are significantly below the strongest winds estimated to occur on this moon for a large range of particle sizes. On Pluto, the surface pressure and temperature conditions measured by New Horizons 2015 \cite{Sternetal15} result in a threshold that is at the limit of aeolian sediment transport for current estimates of the strongest winds (Methods). 
However, seasonal changes induced by its eccentric orbit and extreme axial tilt ($\approx120^\circ$) are likely to increase Pluto's surface pressure via sublimation of nitrogen from about $1$Pa in 2015 to at least $10$Pa in the next 10-20 years \cite{Sternetal16}. By then Pluto's surface conditions are more than sufficient to sustain aeolian sediment transport (Fig.~\ref{Planetpredictions}). 

The prediction for comet 67P has to be interpreted with caution because gravity on 67P decreases rapidly with the distance to its core, whereas our model assumed a constant gravity acceleration. Therefore, steady sediment transport is only possible as long as the characteristic travelling height $H$ and length $L\gg H$ of fast transported particles (Methods) are much smaller than the size of the comet ($\approx2$km). This requirement is only fulfilled for a small range of particle sizes (left of black dashed line in Fig.~3). 

In summary, the predicted sediment transport cessation thresholds suggest active aeolian sediment transport on Triton, Pluto and comet 67P, which could explain the origin of reported surface features. Similarly to what happens on Mars \cite{Kok10a}, the constant ejection of particles from the surface due to particle-bed impacts allows sediment transport to be sustained under relatively weak winds. In fact, estimated winds are usually far too weak to eject surface particles directly, and a few particles must initially enter the air flow by a mechanism distinct from wind erosion to trigger sediment transport. On comet 67P, Rosetta's observations evidenced that surface particles are continuously ejected and subsequently fall back to the surface (``airfall'' \cite{Thomasetal15b}). On Triton, Voyager~2 observed two geyser-like eruptions ejecting large amounts of sediment during its brief flyby \cite{Hansenetal90,SaganChyba90}. And even in the absence of such predominant particle-supplying mechanisms, micrometeoroid impacts eject surface particles at a regular basis, especially in low-atmospheric-density and low-gravity environments, where melting and evaporation of meteorites do not pose much of an issue. We conservatively estimate that about two impacts of micrometeoroids with size $0.1\mathrm{mm}<D<1\mathrm{mm}$ from the Kuiper belt --each such impact nearly guarantees the ejection of surface particles-- occur on Pluto's surface per km$^2$ per year (Methods). This is a large rate considering that the time scale associated with the alteration of Pluto's youngest surface regions (e.g., Sputnik Planum) in the absence of aeolian and fluvial processes is probably of the order of 1-10 million years \cite{Trilling16}.

The recent discovery of the ``Snakeskin'' terrain on Pluto (Fig.~\ref{Fluid_sheared_surfaces}d), which shows ridges, ripples, and perhaps even barchan dunes, led to the speculation that Pluto might have had a much thicker atmosphere millions of years ago \cite{Hand15}. Based on previous estimates of transport threshold, it was argued that only then it could have been possible for the wind to shape these bedforms. Here, we have shown that present-orbit conditions on Pluto are very likely sufficient to alter its surface via aeolian sediment transport.

%% Put the bibliography here, most people will use BiBTeX in
%% which case the environment below should be replaced with
%% the \bibliography{} command.

%\bibliography{model}

%\bibliographystyle{naturemag}
\begin{thebibliography}{10}
\expandafter\ifx\csname url\endcsname\relax
  \def\url#1{\texttt{#1}}\fi
\expandafter\ifx\csname urlprefix\endcsname\relax\def\urlprefix{URL }\fi
\providecommand{\bibinfo}[2]{#2}
\providecommand{\eprint}[2][]{\url{#2}}

\bibitem{Hansenetal90}
\bibinfo{author}{Hansen, C.~J.}, \bibinfo{author}{McEwen, A.~S.},
  \bibinfo{author}{Ingersoll, A.~P.} \& \bibinfo{author}{Terrile, R.~J.}
\newblock \bibinfo{title}{Surface and airborne evidence for plumes and winds on
  triton}.
\newblock \emph{\bibinfo{journal}{Science}} \textbf{\bibinfo{volume}{250}},
  \bibinfo{pages}{421--424} (\bibinfo{year}{1990}).

\bibitem{SaganChyba90}
\bibinfo{author}{Sagan, C.} \& \bibinfo{author}{Chyba, C.}
\newblock \bibinfo{title}{Triton's streaks as windblown dust}.
\newblock \emph{\bibinfo{journal}{Nature}} \textbf{\bibinfo{volume}{346}},
  \bibinfo{pages}{546--548} (\bibinfo{year}{1990}).

\bibitem{Sternetal15}
\bibinfo{author}{Stern, S.~A.} \emph{et~al.}
\newblock \bibinfo{title}{The pluto system: Initial results from its
  exploration by new horizons}.
\newblock \emph{\bibinfo{journal}{Science}} \textbf{\bibinfo{volume}{350}},
  \bibinfo{pages}{aad1815} (\bibinfo{year}{2015}).

\bibitem{Thomasetal15a}
\bibinfo{author}{Thomas, N.} \emph{et~al.}
\newblock \bibinfo{title}{The morphological diversity of comet
  67p/churyumov-gerasimenko}.
\newblock \emph{\bibinfo{journal}{Science}} \textbf{\bibinfo{volume}{347}},
  \bibinfo{pages}{aaa0440} (\bibinfo{year}{2015}).

\bibitem{Thomasetal15b}
\bibinfo{author}{Thomas, N.} \emph{et~al.}
\newblock \bibinfo{title}{Redistribution of particles across the nucleus of
  comet 67p/churyumov-gerasimenko}.
\newblock \emph{\bibinfo{journal}{Astronomy \& Astrophysics}}
  \textbf{\bibinfo{volume}{583}}, \bibinfo{pages}{A17} (\bibinfo{year}{2015}).

\bibitem{Mooreetal15}
\bibinfo{author}{Moore, J.~M.} \emph{et~al.}
\newblock \bibinfo{title}{Geology before pluto: Pre-encounter considerations}.
\newblock \emph{\bibinfo{journal}{Icarus}} \textbf{\bibinfo{volume}{246}},
  \bibinfo{pages}{65--81} (\bibinfo{year}{2015}).

\bibitem{Gladstoneetal16}
\bibinfo{author}{Gladstone, G.~R.} \emph{et~al.}
\newblock \bibinfo{title}{The atmosphere of pluto as observed by new horizons}.
\newblock \emph{\bibinfo{journal}{Science}} \textbf{\bibinfo{volume}{351}},
  \bibinfo{pages}{aad8866} (\bibinfo{year}{2016}).

\bibitem{Kok10a}
\bibinfo{author}{Kok, J.~F.}
\newblock \bibinfo{title}{Difference in the wind speeds required for initiation
  versus continuation of sand transport on mars: Implications for dunes and
  dust storms}.
\newblock \emph{\bibinfo{journal}{Physical Review Letters}}
  \textbf{\bibinfo{volume}{104}}, \bibinfo{pages}{074502}
  (\bibinfo{year}{2010}).

\bibitem{ClaudinAndreotti06}
\bibinfo{author}{Claudin, P.} \& \bibinfo{author}{Andreotti, B.}
\newblock \bibinfo{title}{A scaling law for aeolian dunes on mars, venus,
  earth, and for subaqueous ripples}.
\newblock \emph{\bibinfo{journal}{Earth and Planetary Science Letters}}
  \textbf{\bibinfo{volume}{252}}, \bibinfo{pages}{30--44}
  (\bibinfo{year}{2006}).

\bibitem{Kok10b}
\bibinfo{author}{Kok, J.~F.}
\newblock \bibinfo{title}{An improved parameterization of wind blown sand flux
  on mars that includes the effect of hysteresis}.
\newblock \emph{\bibinfo{journal}{Geophysical Research Letters}}
  \textbf{\bibinfo{volume}{37}}, \bibinfo{pages}{L12202}
  (\bibinfo{year}{2010}).

\bibitem{Paehtzetal12}
\bibinfo{author}{P\"ahtz, T.}, \bibinfo{author}{Kok, J.~F.} \&
  \bibinfo{author}{Herrmann, H.~J.}
\newblock \bibinfo{title}{The apparent roughness of a sand surface blown by
  wind from an analytical model of saltation}.
\newblock \emph{\bibinfo{journal}{New Journal of Physics}}
  \textbf{\bibinfo{volume}{14}}, \bibinfo{pages}{043035}
  (\bibinfo{year}{2012}).

\bibitem{Duranetal12}
\bibinfo{author}{Dur\'an, O.}, \bibinfo{author}{Andreotti, B.} \&
  \bibinfo{author}{Claudin, P.}
\newblock \bibinfo{title}{Numerical simulation of turbulent sediment transport,
  from bed load to saltation}.
\newblock \emph{\bibinfo{journal}{Physics of Fluids}}
  \textbf{\bibinfo{volume}{24}}, \bibinfo{pages}{103306}
  (\bibinfo{year}{2012}).

\bibitem{PaehtzDuran16a}
\bibinfo{author}{P\"ahtz, T.} \& \bibinfo{author}{Dur\'an, O.}
\newblock \bibinfo{title}{Fluid forces or impacts, what governs the entrainment
  of soil particles in sediment transport mediated by a newtonian fluid?}
\newblock \emph{\bibinfo{journal}{under review in Physical Review Letters}}
  (\bibinfo{year}{2016}).

\bibitem{Clarketal15}
\bibinfo{author}{Clark, A.~H.}, \bibinfo{author}{Shattuck, M.~D.},
  \bibinfo{author}{Ouellette, N.~T.} \& \bibinfo{author}{O'Hern, C.~S.}
\newblock \bibinfo{title}{Onset and cessation of motion in hydrodynamically
  sheared granular beds}.
\newblock \emph{\bibinfo{journal}{Physical Review E}}
  \textbf{\bibinfo{volume}{92}}, \bibinfo{pages}{042202}
  (\bibinfo{year}{2015}).

\bibitem{Maurinetal15}
\bibinfo{author}{Maurin, R.}, \bibinfo{author}{Chauchat, J.},
  \bibinfo{author}{Chareyre, B.} \& \bibinfo{author}{Frey, P.}
\newblock \bibinfo{title}{A minimal coupled fluid-discrete element model for
  bedload transport}.
\newblock \emph{\bibinfo{journal}{Physics of Fluids}}
  \textbf{\bibinfo{volume}{27}}, \bibinfo{pages}{113302}
  (\bibinfo{year}{2015}).

\bibitem{Koketal12}
\bibinfo{author}{Kok, J.~F.}, \bibinfo{author}{Parteli, E. J.~R.},
  \bibinfo{author}{Michaels, T.~I.} \& \bibinfo{author}{Karam, D.~B.}
\newblock \bibinfo{title}{The physics of wind-blown sand and dust}.
\newblock \emph{\bibinfo{journal}{Reports on Progress in Physics}}
  \textbf{\bibinfo{volume}{75}}, \bibinfo{pages}{106901}
  (\bibinfo{year}{2012}).

\bibitem{Bagnold56}
\bibinfo{author}{Bagnold, R.~A.}
\newblock \bibinfo{title}{The flow of cohesionless grains in fluid}.
\newblock \emph{\bibinfo{journal}{Philosophical Transactions of the Royal
  Society London A}} \textbf{\bibinfo{volume}{249}}, \bibinfo{pages}{235--297}
  (\bibinfo{year}{1956}).

\bibitem{Camenen07}
\bibinfo{author}{Camenen, B.}
\newblock \bibinfo{title}{Simple and general formula for the settling velocity
  of particles}.
\newblock \emph{\bibinfo{journal}{Journal of Hydraulic Engineering}}
  \textbf{\bibinfo{volume}{133}}, \bibinfo{pages}{229--233}
  (\bibinfo{year}{2007}).

\bibitem{Croweetal12}
\bibinfo{author}{Crowe, C.~T.}, \bibinfo{author}{Schwarzkopf, J.~D.},
  \bibinfo{author}{Sommerfeld, M.} \& \bibinfo{author}{Tsuji, Y.}
\newblock \emph{\bibinfo{title}{Multiphase Flows with Droplets and Particles}}
  (\bibinfo{publisher}{Taylor \& Francis Group, Boca Raton},
  \bibinfo{year}{2012}).

\bibitem{Buffingtonmontgomery97}
\bibinfo{author}{Buffington, J.~M.} \& \bibinfo{author}{Montgomery, D.~R.}
\newblock \bibinfo{title}{A systematic analysis of eight decades of incipient
  motion studies, with special reference to gravel-bedded rivers}.
\newblock \emph{\bibinfo{journal}{Water Resources Research}}
  \textbf{\bibinfo{volume}{33}}, \bibinfo{pages}{1993--2029}
  (\bibinfo{year}{1997}).

\bibitem{Ouremietal07}
\bibinfo{author}{Ouriemi, M.}, \bibinfo{author}{Aussillous, P.},
  \bibinfo{author}{Medale, M.}, \bibinfo{author}{Peysson, Y.} \&
  \bibinfo{author}{Guazzelli, E.}
\newblock \bibinfo{title}{Determination of the critical shields number for
  particle erosion in laminar flow}.
\newblock \emph{\bibinfo{journal}{Physics of Fluids}}
  \textbf{\bibinfo{volume}{19}}, \bibinfo{pages}{061706}
  (\bibinfo{year}{2007}).

\bibitem{Bagnold37}
\bibinfo{author}{Bagnold, R.~A.}
\newblock \bibinfo{title}{The transport of sand by wind}.
\newblock \emph{\bibinfo{journal}{The Geographical Journal}}
  \textbf{\bibinfo{volume}{89}}, \bibinfo{pages}{409--438}
  (\bibinfo{year}{1937}).

\bibitem{Chepil45}
\bibinfo{author}{Chepil, W.~S.}
\newblock \bibinfo{title}{Dynamics of wind erosion: Ii. initiation of soil
  movement}.
\newblock \emph{\bibinfo{journal}{Soil Science}} \textbf{\bibinfo{volume}{60}},
  \bibinfo{pages}{397--411} (\bibinfo{year}{1945}).

\bibitem{Creysselsetal09}
\bibinfo{author}{Creyssels, M.} \emph{et~al.}
\newblock \bibinfo{title}{Saltating particles in a turbulent boundary layer:
  experiment and theory}.
\newblock \emph{\bibinfo{journal}{Journal of Fluid Mechanics}}
  \textbf{\bibinfo{volume}{625}}, \bibinfo{pages}{47--74}
  (\bibinfo{year}{2009}).

\bibitem{Hoetal11}
\bibinfo{author}{Ho, T.~D.}, \bibinfo{author}{Valance, A.},
  \bibinfo{author}{Dupont, P.} \& \bibinfo{author}{{Ould El Moctar}, A.}
\newblock \bibinfo{title}{Scaling laws in aeolian sand transport}.
\newblock \emph{\bibinfo{journal}{Physical Review Letters}}
  \textbf{\bibinfo{volume}{106}}, \bibinfo{pages}{094501}
  (\bibinfo{year}{2011}).

\bibitem{Sternetal16}
\bibinfo{author}{Stern, S.~A.} \emph{et~al.}
\newblock \bibinfo{title}{Evidence that pluto's long-term polar axis shifts
  drive sharp changes in the planet’s atmospheric pressure over time}.
\newblock \emph{\bibinfo{journal}{under review in Icarus}}
  (\bibinfo{year}{2016}).
\newblock \bibinfo{note}{Video (`Press Briefing - New Horizons: Peering into
  Pluto’s Past'): http://livestream.com/viewnow/LPSC2016}.

\bibitem{Trilling16}
\bibinfo{author}{Trilling, D.~E.}
\newblock \bibinfo{title}{The surface age of sputnik planum, pluto, must be
  less than 10 million years}.
\newblock \emph{\bibinfo{journal}{PLOS ONE}} \textbf{\bibinfo{volume}{11}},
  \bibinfo{pages}{e0147386} (\bibinfo{year}{2016}).

\bibitem{Hand15}
\bibinfo{author}{Hand, E.}
\newblock \bibinfo{title}{Late harvest from pluto reveals a complex world}.
\newblock \emph{\bibinfo{journal}{Science}} \textbf{\bibinfo{volume}{350}},
  \bibinfo{pages}{260--261} (\bibinfo{year}{2015}).

\bibitem{Burretal15}
\bibinfo{author}{Burr, D.~M.} \emph{et~al.}
\newblock \bibinfo{title}{Higher-than-predicted saltation threshold wind speeds
  on titan}.
\newblock \emph{\bibinfo{journal}{Nature}} \textbf{\bibinfo{volume}{517}},
  \bibinfo{pages}{60--63} (\bibinfo{year}{2015}).

\bibitem{Ayoubetal14}
\bibinfo{author}{Ayoub, F.} \emph{et~al.}
\newblock \bibinfo{title}{Threshold for sand mobility on mars calibrated from
  seasonal variations of sand flux}.
\newblock \emph{\bibinfo{journal}{Nature Communication}}
  \textbf{\bibinfo{volume}{5}}, \bibinfo{pages}{5096} (\bibinfo{year}{2014}).

\bibitem{Tokano10}
\bibinfo{author}{Tokano, T.}
\newblock \bibinfo{title}{Relevance of fast westerlies at equinox for the
  eastward elongation of titan's dunes}.
\newblock \emph{\bibinfo{journal}{Aeolian Research}}
  \textbf{\bibinfo{volume}{2}}, \bibinfo{pages}{113--127}
  (\bibinfo{year}{2010}).

\bibitem{Fentonetal05}
\bibinfo{author}{Fenton, L.~K.}, \bibinfo{author}{Toigo, A.~D.} \&
  \bibinfo{author}{Richardson, M.~I.}
\newblock \bibinfo{title}{Aeolian processes in proctor crater on mars:
  Mesoscale modeling of dune-forming winds}.
\newblock \emph{\bibinfo{journal}{Journal of Geophysical Research}}
  \textbf{\bibinfo{volume}{110}}, \bibinfo{pages}{E06005}
  (\bibinfo{year}{2005}).

\bibitem{Bagnold38}
\bibinfo{author}{Bagnold, R.~A.}
\newblock \bibinfo{title}{The measurement of sand storms}.
\newblock \emph{\bibinfo{journal}{Proceedings of the Royal Society London
  Series A}} \textbf{\bibinfo{volume}{167}}, \bibinfo{pages}{282--291}
  (\bibinfo{year}{1938}).

\bibitem{GreeleyArvidson90}
\bibinfo{author}{Greeley, R.} \& \bibinfo{author}{Arvidson, R.~E.}
\newblock \bibinfo{title}{Aeolian processes on venus}.
\newblock \emph{\bibinfo{journal}{Earth, Moon, and Planets}}
  \textbf{\bibinfo{volume}{50}}, \bibinfo{pages}{127--157}
  (\bibinfo{year}{1990}).

\bibitem{Greenstreetetal15}
\bibinfo{author}{Greenstreet, S.}, \bibinfo{author}{Gladman, B.} \&
  \bibinfo{author}{McKinnon, W.~B.}
\newblock \bibinfo{title}{Impact and cratering rates onto pluto}.
\newblock \emph{\bibinfo{journal}{Icarus}} \textbf{\bibinfo{volume}{258}},
  \bibinfo{pages}{267--288} (\bibinfo{year}{2015}).

\bibitem{Zahnleetal03}
\bibinfo{author}{Zahnle, K.}, \bibinfo{author}{Schenk, P.},
  \bibinfo{author}{Levison, H.} \& \bibinfo{author}{Donesc, L.}
\newblock \bibinfo{title}{Cratering rates in the outer solar system}.
\newblock \emph{\bibinfo{journal}{Icarus}} \textbf{\bibinfo{volume}{163}},
  \bibinfo{pages}{263--289} (\bibinfo{year}{2003}).

\bibitem{FlynnMcKay90}
\bibinfo{author}{Flynn, G.~J.} \& \bibinfo{author}{McKay, D.~S.}
\newblock \bibinfo{title}{An assessment of the meteoritic contribution to the
  martian soil}.
\newblock \emph{\bibinfo{journal}{Journal of Geophysical Research}}
  \textbf{\bibinfo{volume}{95}}, \bibinfo{pages}{14497--14509}
  (\bibinfo{year}{1990}).

\bibitem{McKinnonetal08}
\bibinfo{author}{McKinnon, W.~B.}, \bibinfo{author}{Prialnik, D.},
  \bibinfo{author}{Stern, S.~A.} \& \bibinfo{author}{Coradini, A.}
\newblock \bibinfo{title}{Structure and evolution of kuiper belt objects and
  dwarf planets}.
\newblock In \bibinfo{editor}{Barucci, M.~A.}, \bibinfo{editor}{Boehnhardt,
  H.}, \bibinfo{editor}{Cruikshank, D.~P.} \& \bibinfo{editor}{Morbidelli, A.}
  (eds.) \emph{\bibinfo{booktitle}{The Solar System Beyond Neptune}},
  \bibinfo{pages}{213--241} (\bibinfo{publisher}{University of Arizona Press,
  Tucson}, \bibinfo{year}{2008}).

\bibitem{Paehtzetal15a}
\bibinfo{author}{P\"ahtz, T.}, \bibinfo{author}{Dur\'an, O.},
  \bibinfo{author}{Ho, T.-D.}, \bibinfo{author}{Valance, A.} \&
  \bibinfo{author}{Kok, J.~F.}
\newblock \bibinfo{title}{The fluctuation energy balance in non-suspended
  fluid-mediated particle transport}.
\newblock \emph{\bibinfo{journal}{Physics of Fluids}}
  \textbf{\bibinfo{volume}{27}}, \bibinfo{pages}{013303}
  (\bibinfo{year}{2015}).

\bibitem{Carneiroetal11}
\bibinfo{author}{Carneiro, M.~V.}, \bibinfo{author}{P\"ahtz, T.} \&
  \bibinfo{author}{Herrmann, H.~J.}
\newblock \bibinfo{title}{Jump at the onset of saltation}.
\newblock \emph{\bibinfo{journal}{Physical Review Letters}}
  \textbf{\bibinfo{volume}{107}}, \bibinfo{pages}{098001}
  (\bibinfo{year}{2011}).

\bibitem{GuoJulien07}
\bibinfo{author}{Guo, J.} \& \bibinfo{author}{Julien, P.~Y.}
\newblock \bibinfo{title}{Buffer law and transitional roughness effect in
  turbulent open-channel flows}.
\newblock In \emph{\bibinfo{booktitle}{The Fifth International Symposium on
  Environmental Hydraulics (ISEH V), Tempe, Arizona}}, vol.
  \bibinfo{volume}{Paper 5}, \bibinfo{pages}{1--6}
  (\bibinfo{publisher}{University of Nebraska - Lincoln},
  \bibinfo{year}{2007}).

\bibitem{IshiiHibiki11}
\bibinfo{author}{Ishii, M.} \& \bibinfo{author}{Hibiki, T.}
\newblock \emph{\bibinfo{title}{of Two-Phase Flow Thermo-Fluid Dynamics, 323
  pp.}} (\bibinfo{publisher}{Springer, New York}, \bibinfo{year}{2011}).

\end{thebibliography}


%% Here is the endmatter stuff: Supplementary Info, etc.
%% Use \item's to separate, default label is "Acknowledgements"

\begin{addendum}
 \item We thank Pan Jia, Bruno Andreotti, and Philippe Claudin for providing us the surface pressure and temperature values of comet 67P that they obtained from their coma model. We acknowledge support from grants National Natural Science Foundation of China (Nos. 1151101041 and 41376095) and Natural Science Foundation of Zhejiang Province (No. LR16E090001).
 \item[Competing Interests] The authors declare that they have no
competing financial interests.
 \item[Correspondence] Correspondence and requests for materials
should be addressed to T. P\"ahtz~(email: 0012136@zju.edu.cn).
\end{addendum}
\newpage

%%
%% TABLES
%%
%% If there are any tables, put them here.
%%
\setcounter{equation}{0}
\setcounter{figure}{0}
\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thefigure}{(Extended Data) \arabic{figure}}

\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Paehtz_EDfig1.pdf}
 \end{center}
 \caption{ {\bf Model comparison}. Cessation thresholds predicted by this and previous models based on the particle-bed-impact entrainment mechanism~\cite{ClaudinAndreotti06,Kok10b,Paehtzetal12} (solid lines). A recent semi-empirical prediction of the initiation threshold \cite{Burretal15} is also included for reference (dashed line). Models~\cite{ClaudinAndreotti06,Kok10b,Paehtzetal12} have been slightly modified to account for the buoyancy and added-mass force. Previous models strongly overestimate measurements in water even though turbulent subaqueous sediment transport is sustained through impacts \cite{Clarketal15,PaehtzDuran16a}.}
\end{figure}

\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Paehtz_EDfig2.pdf}
 \end{center}
 \caption{{\bf Validation of Eq.~(1).} Cohesionless friction coefficient $\mu^o_b=\mu_b(0)$ versus the Shields parameter $\Theta$. Symbols correspond to data from transport simulations~\cite{Duranetal12} for varying $s$, $\mathrm{Re}_p$, and $\Theta$. The solid line corresponds to the constant value $\mu^o_b=0.67$ assumed in the analytical model.}
\end{figure}

\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Paehtz_EDfig3.pdf}
 \end{center}
 \caption{{\bf Validation of Eq.~(3).} (a) $V_x/\sqrt{\Theta^{\mathrm{ex}}_t}$ versus $Z/Z_\nu$. (b) $V_x/\sqrt{\Theta^{\mathrm{ex}}_t}$ versus $U_x/\sqrt{\Theta^{\mathrm{ex}}_t}$. Symbols correspond to data from transport simulations~\cite{Duranetal12} for varying $s$, $\mathrm{Re}_p$, and $\Theta$ close to the threshold. The solid line corresponds to perfect agreement.}
\end{figure}

\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Paehtz_EDfig4.pdf}
 \end{center}
 \caption{{\bf Validation of Eq.~(4).} Dimensionless mean transport layer height $Z$ versus $Z_c+sV_z^2$. Symbols correspond to data from transport simulations~\cite{Duranetal12} for varying $s$, $\mathrm{Re}_p$, and $\Theta$. The solid line corresponds to perfect agreement.}
\end{figure}

\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Paehtz_EDfig5.pdf}
 \end{center}
 \caption{{\bf Validation of Eq.~(5).} (a) $\overline{|v_z|}/\sqrt{(s-1)gd}$ versus $\overline{\mu}^{-1}\sqrt{V_x|_{\Theta\approx\Theta^{\mathrm{ex}}_t}V_x}$ and (b) $V_z \equiv \sqrt{\overline{v^2_z}}/\sqrt{(s-1)gd}$ versus $\mu_b^{-1}\sqrt{V_x|_{\Theta\approx\Theta^{\mathrm{ex}}_t}V_x}$. Symbols correspond to data from transport simulations~\cite{Duranetal12} for varying $s$, $\mathrm{Re}_p$, and $\Theta$. The solid line corresponds to (a) Eq.~(\ref{beforeclosure}) and (b) Eq.~(5) with $V_x\rightarrow\sqrt{V_x|_{\Theta\approx\Theta^{\mathrm{ex}}_t}V_x}$ and $\alpha=0.18$.}
\end{figure}

\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Paehtz_EDfig6.pdf}
 \end{center}
 \caption{{\bf Test of the final threshold prediction with cohesionless numerical particle-scale simulations.} Predicted versus simulated values of $\Theta^{\mathrm{ex}}_t$ for varying $s$ and $\mathrm{Re}_p$.}
\end{figure}

\begin{figure}[htb]
 \begin{center}
  \includegraphics[width=1.0\columnwidth]{Paehtz_EDfig7.pdf}
 \end{center}
 \caption{{\bf Validation of Eqs.~(\ref{Visc}-\ref{Turb}).} Analytical predictions of the threshold Shields parameter $\Theta^{\mathrm{ex}}_t$ for non-rarefied ($\mathrm{Kn}=0$), non-cohesive ($C=0$) conditions as a function of the impact number $\mathrm{Im}$ for different values of $s$. Solid lines correspond to predictions from our model, whereas dashed lines correspond to predictions from the model of Ref.~\cite{ClaudinAndreotti06} (slightly modified to account for the buoyancy and added-mass force).}
\end{figure}

\newpage
\section*{Methods}
\subsection{Environmental conditions}
\subsubsection*{Water}
We assume that particles under water consist of quartz ($\rho_p=2650$kg/m$^3$, $E=7\times10^{10}$Pa, $\gamma=3$J/m$^2$), and that water has a density $\rho_f=1000$kg/m$^3$, corresponding to $s=2.65$, $\nu=1\times10^{-6}$m$^2$/s, and speed of sound $c\approx1500$m/s. The gravitational constant is taken as $g=9.81$m/s$^2$.

\subsubsection*{Venus}
We assume that Venus' particles consist of basalt ($\rho_p=3000$kg/m$^3$, $E=7\times10^{10}$Pa, $\gamma=3$J/m$^2$) and Venus' atmosphere to $96.5\%$ of CO$_2$ and $3.5\%$ of N$_2$ with surface pressure $P=9.2\times10^6$Pa and temperature $T=737$K. This corresponds to $s\approx45$, $\nu\approx4.9\times10^{-7}$m$^2$/s (cf. \cite{Burretal15}), and $c=440$m/s. The gravitational constant is taken as $g=8.87$m/s$^2$.

\subsubsection*{Titan}
We assume that Titan's particles consist of water ice and/or organics ($\rho_p=1000$kg/m$^3$, $E=9\times10^9$Pa, $\gamma=0.11$J/m$^2$), and that Titan's atmosphere has a density $\rho_f=5.3$kg/m$^3$, corresponding to $s=190$, $\nu=1.2\times10^{-6}$m$^2$/s \cite{Burretal15}, and $c\approx200$m/s. The gravitational constant is taken as $g=1.35$m/s$^2$.

\subsubsection*{Earth}
We assume that Earth's particles consist of quartz or clay loam as in the experiments of Refs.~\cite{Bagnold37,Chepil45} ($\rho_p=1650-2650$kg/m$^3$, $E=70\times10^9$Pa, $\gamma=3$J/m$^2$), and that Earth's atmosphere has a density $\rho_f=1.2$kg/m$^3$, corresponding to $s=1350-2210$, $\nu=1.6\times10^{-5}$m$^2$/s, and $c\approx340$m/s. Since most particles in the experiments are in the upper end of the density range, we choose $s=2000$ ($\rho_p=2400$kg/m$^3$). Note that the choice of $s$ does not much affect the agreement with the experiments. The gravitational constant is taken as $g=9.81$m/s$^2$.

\subsubsection*{Mars}
We assume that Mars' particles consist of basalt ($\rho_p=3000$kg/m$^3$, $E=7\times10^{10}$Pa, $\gamma=3$J/m$^2$) and Mars' atmosphere to $96\%$ of CO$_2$ and $2\%$ of N$_2$ and $2\%$ of Ar with $P=520$Pa and $T=230$K \cite{Ayoubetal14}. This corresponds to $s\approx2.5\times10^5$, $\nu\approx9.9\times10^{-4}$m$^2$/s, and $c\approx250$m/s. The gravitational constant is taken as $g=3.71$m/s$^2$.

\subsubsection*{Triton}
We assume that Triton's particles consist of methane ice ($\rho_p=500$kg/m$^3$, $E=9\times10^9$Pa, $\gamma=0.11$J/m$^2$) and Triton's atmosphere of N$_2$ with $P=1.6$Pa and $T=38$K \cite{SaganChyba90}. This corresponds to $s\approx3.5\times10^6$, $\nu\approx0.016$m$^2$/s, and $c\approx130$m/s. The gravitational constant is taken as $g=0.78$m/s$^2$.

\subsubsection*{Pluto}
We assume that Pluto's particles consist of methane ice ($\rho_p=500$kg/m$^3$, $E=9\times10^9$Pa, $\gamma=0.11$J/m$^2$) and Pluto's atmosphere of N$_2$ with $P=1$Pa and $T=38$K \cite{Sternetal15}. This corresponds to $s\approx5.6\times10^6$, $\nu\approx0.025$m$^2$/s, and $c\approx130$m/s. When Pluto reaches its largest-pressure conditions in the current orbit, we assume $P=10$Pa and T=$43$K. This corresponds to $s\approx6.5\times10^5$, $\nu\approx0.0033$m$^2$/s, and $c\approx130$m/s. The gravitational constant is taken as $g=0.62$m/s$^2$.

\subsubsection*{Comet 67P/Churyumov-Gerasimenko}
Philippe Claudin and co-authors have developed a model for the coma of 67P (consisting of H$_2$O). We use the values they obtained from their model at perihelion in the morning ($P=0.15$Pa, $T=200K$, private correspondence with Philippe Claudin). Assuming $\rho_p\approx1000kg$/m$^3$, $E=9\times10^9$Pa, and $\gamma=0.11$J/m$^2$ (organics), this corresponds to $s\approx6.2\times10^{8}$, $\nu\approx4.9$m$^2$/s, and $c\approx360$m/s. The gravitational constant is taken as $g=1.55\times10^{-4}$m/s$^2$ \cite{Thomasetal15b}.

\subsubsection*{Strongest winds on planetary bodies}
Ref.~\cite{SaganChyba90} estimated a maximal wind shear velocity $u_\ast=0.75$m/s on Triton from a maximal geostrophic wind speed $u_\infty=16.7$m/s via $u_\ast^2=C_d^\mathrm{bound}u_\infty^2$, where $C_d^\mathrm{bound}\approx0.002$ is the boundary layer drag coefficient. Using the same method, we estimate a maximal wind shear velocity of about $0.45$m/s on Pluto from a geostrophic wind speed of $10$m/s, recently estimated as an upper bound by Ref.~\cite{Gladstoneetal16}. On Titan, general circulation models suggest free stream velocities of up to $1.5$m/s \cite{Tokano10}, corresponding to wind shear velocities of up to $0.067$m/s (same method). On Mars, general circulation models suggest wind shear velocities of up to $1.2$m/s \cite{Fentonetal05}. For Earth, we assume a maximal value of $0.82$m/s, which was measured during a strong sand storm \cite{Bagnold38}. For Venus, we estimate a maximal wind shear velocity $u_\ast=0.12$m/s from measured maximal $1$m winds of about $2$m/s \cite{GreeleyArvidson90} using a logarithmic wind profile $u_\ast/\kappa\ln z/z_o$ assuming a surface roughness $z_o=1$mm. Moreover, using numerical simulations, Ref.~\cite{Thomasetal15b} estimated maximal wind shear velocities of about $500$m/s on comet 67P at perihelion in the region where the ripples are located (the same value is obtained for perihelion in the morning from the coma model by Philippe Claudin and co-workers - private correspondence with Philippe Claudin). The reason for such extreme winds is the irregular shape of the nucleus that, combined with the lateral expansion toward the nightside, partially funnels the gas over the ripple area.

\subsubsection*{Cosmic impact rate on Pluto}
Ref.~\cite{Greenstreetetal15} inferred a total rate of $4.8\times10^{-11}$ impacts on Pluto of Kuiper belt objects with $D>108$km, corresponding to a rate of $1.9\times10^{-16}$ impacts on Pluto per object per year (see their Table~1). Since older studies reported larger values ($5{-}7\times10^{-16}$) when smaller particles are taken into account \cite{Zahnleetal03}, we take the value of Ref.~\cite{Greenstreetetal15} as a conservative estimate. Ref.~\cite{Greenstreetetal15} also reviewed studies on the size distribution of Kuiper belt objects with $D\lessapprox100$km, very roughly suggesting $\d N\propto D^{-3}$. Hence, the rate of encounters of Pluto with objects with $D_-<D<D_+$ per year is estimated as $4.8\times10^{-11}[(D_-/108\mathrm{km})^{-2}-(D_+/108\mathrm{km})^{-2}]$. For $D_-=0.1$mm and $D_+=1$mm, this corresponds to $5.5\times10^7$ encounters per year. A large fraction of rocky objects with this size range already survive the passage through Mars' atmosphere without being melted due to its large atmospheric scale height and small escape velocity (i.e., the relevant entry velocity scale for meteoroids captured by the gravity of the target) compared with Earth \cite{FlynnMcKay90}. On Pluto, the former is even much larger and the latter even much smaller. Hence, we can safely assume that all rocky micrometeorites survive the passage through Pluto's atmosphere and only the icy ones are being melted. Assuming an ice-to-rock ratio of $30/70$ for Kuiper belt objects \cite{McKinnonetal08} and a homogeneous distribution of impact locations on Pluto's surface, the encounter rate of micrometeorites with $0.1\mathrm{mm}<D<1\mathrm{mm}$ estimated above corresponds to about two impacts per km$^2$ per year.

\subsection{Model}

\subsubsection*{Definition of local average over particles}
The particle-mass-weighted average of a quantity $A$ over all particles with the same elevation $z$ is defined as \cite{Paehtzetal15a}
\begin{eqnarray}
 \langle A\rangle=\frac{1}{\rho}\overline{\sum_nm^nA^n\delta(\mathbf{x}-\mathbf{x}^n)}^E,
\end{eqnarray}
where $\mathbf{x}^{(n)}$ is the location (of particle $n$), $\rho=\overline{\sum_nm^n\delta(\mathbf{x}-\mathbf{x}^n)}^E$ the particle mass density, $\delta$ the delta distribution, and the overbar with superscript `E' denotes the ensemble average.

\subsubsection*{Definition of transport layer average}

The average over the transport layer of a quantity $A$ is represented by an overbar and defined as
\begin{eqnarray}
 \overline{A}=\int\limits_{z_s}^\infty\rho\langle A\rangle\d z/\int\limits_{z_s}^\infty\rho\d z,
\end{eqnarray}
where $z_s$ is the bed surface elevation.

\subsubsection*{Precise definition of $\Theta^{\mathrm{ex}}_t$}
We define $\Theta^{\mathrm{ex}}_t$ from extrapolating $M=Q^2/\int_{-\infty}^\infty\rho\langle v_x\rangle^2\d z$ to zero via $M\propto\Theta-\Theta^{\mathrm{ex}}_t$, where $Q=\int_{-\infty}^\infty\rho\langle v_x\rangle\d z$ is the sediment transport rate. This method is almost identical to the method described by Ref.~\cite{Duranetal12} and yields nearly the same values. For aeolian conditions, extrapolating $Q$ to zero via $Q\propto\Theta-\Theta^{\mathrm{ex}}_t$ also yields almost the same values in our simulations. We use the latter method to obtain $\Theta^{\mathrm{ex}}_t$ from the experiments by Refs.~\cite{Creysselsetal09,Hoetal11}. Note that $\Theta^{\mathrm{ex}}_t$ is smaller than the threshold $\Theta_t$ below which sediment transport ceases. That is, the sediment transport rate $Q$ (nearly) vanishes when $\Theta<\Theta_t$ and jumps to a finite value at $\Theta_t$ \cite{Carneiroetal11,Clarketal15}. However, this does not pose an issue as our simulations indicate that $\Theta^{\mathrm{ex}}_t$ and $\Theta_t$ are nearly identical when $\mathrm{Re}_p<10$. On Triton, Pluto, and comet 67P, $\mathrm{Re}_p\ll10$ for the conditions that suggest transport (shaded areas in Fig.~3).

\subsubsection*{Definition of (bed) friction coefficient}

The friction coefficient is defined as
\begin{eqnarray}
 \mu=-\frac{P_{zx}}{P_{zz}},
\end{eqnarray}
where $P_{ij}$ is the particle stress tensor composed of a transport contribution (superscript `t') and a contact contribution (superscript `c'), reading \cite{Paehtzetal15a}
\begin{eqnarray}
 P_{ij}&=&P^t_{ij}+P^c_{ij} \label{Particlestress}, \\
 P^t_{ij}&=&\rho(\langle v_iv_j\rangle-\langle v_i\rangle\langle v_j\rangle), \\
 P^c_{ij}&=&\frac{1}{2}\overline{\sum_{mn}F^{mn}_j(x^m_i-x^n_i)\int\limits_0^1\delta(\mathbf{x}-((\mathbf{x}^m-\mathbf{x}^n)s^\prime+\mathbf{x}^n))\d s^\prime}^E,
\end{eqnarray}
where $\mathbf{F}^{mn}$ is the contact force applied by particle $n$ on particle $m$ and $s^\prime$ a dummy variable. Note that $P^t_{zi}=\rho\langle v_zv_i\rangle$ since $\langle v_z\rangle=0$. The bed friction coefficient is defined as $\mu_b=\mu(z_s)$.

\subsubsection*{Definition of bed surface}

We define the bed surface $z_s$ as the elevation where the granular pressure work is maximal \cite{PaehtzDuran16a}. Granular shear work is defined as $P_{zz}\d\langle v_x\rangle/\d z$.

\subsubsection*{Average drag acceleration (function $f_1$ in Eq.~(1))}

Averaging the horizontal drag acceleration $a^d_x(z)$ over the transport layer leads to,
\begin{equation}
	\frac{\overline{a^d_x}}{\tilde g}= f_1(\mathrm{Re}_p ,U_x-V_x,V_z,\mathrm{Kn}) = \frac{3}{4}C_d^{\rm eff}(U_x-V_x)^2.
\end{equation}
The effective drag coefficient $C_d^{\rm eff}$ includes the effects of rarefaction and is a function of the fluid Reynolds number, relative to the particles, $\mathrm{Re}_r=\mathrm{Re}_p \sqrt{(U_x-V_x)^2+2V_z^2}$ and the relative Mach number $\mathrm{M}_r = \sqrt{(U_x-V_x)^2+2V_z^2}/[c/\sqrt{(s-1)gd}]$ (assuming isotropic particle velocity fluctuations, see below), where $c$ is the speed of sound, and the other quantities are defined in the text. Note that the dependence on the Knudsen number $\mathrm{Kn}$ in Eq.~(1) is parametrized by the ratio $\mathrm{M}_r/\mathrm{Re}_r$. From Ref.~\cite{Croweetal12},
\begin{eqnarray}
 &&C_d^{\rm eff} = 2 C_{\rm fluc}+\left( C_d-2 C_{\rm fluc}\right)\exp\left( -3.07\sqrt{k}\,h_1(\mathrm{Re}_r)\frac{\mathrm{M}_r}{\mathrm{Re}_r}\right)+\frac{h_2(\mathrm{M}_r)}{\sqrt{k}\mathrm{M}_r}\exp\left(-\frac{\mathrm{Re}_r}{2\mathrm{M}_r}\right), \nonumber \\
 &&h_1(\mathrm{Re}_r)=\frac{1+\mathrm{Re}_r(12.278+0.548\mathrm{Re}_r)}{1+11.278\mathrm{Re}_r}, \nonumber \\
 &&h_2(\mathrm{M}_r)=\frac{5.6}{1+\mathrm{M}_r}+1.7, \label{Dragrare}
\end{eqnarray}
where $k$ is the adiabatic exponent, which is usually between about $1.3$ and $1.6$ (we choose $k=1.4$), 
$C_{\rm fluc}$ is the drag enhancement factor due to particle velocity fluctuations (see below), and $C_d$ is the average drag coefficient, which combines viscous and turbulent effects~\cite{Camenen07},
\begin{eqnarray}
 C_d = \left[\left(\frac{24}{\mathrm{Re}_p (U_x-V_x)}\right)^{1/m}+ (C_{\rm fluc}C_d^\infty)^{1/m}\right]^m,
\end{eqnarray}
where $C_d^\infty$ is the particle turbulent drag coefficient and $m$ is the drag exponent ($C_d^\infty=1$ and $m=1.5$ for natural particles and $C_d^\infty=0.5$ and $m=2$ for spheres as in the transport simulations). The factor $C_{\rm fluc}$ is obtained after averaging the horizontal turbulent drag acceleration of a single particle \cite{Camenen07} $a^d_x|_\mathrm{turb}=\frac{3}{4sd}C_d^\infty|\mathbf{u}-\mathbf{v}|(u_x-v_x)$:
\begin{eqnarray}
 \langle a^d_x\rangle |_{\rm turb}=\frac{3}{4sd}C_d^\infty\left[\sqrt{\langle u_x-v_x\rangle^2+\langle v_y^2+v_z^2\rangle}\langle u_x-v_x\rangle+\langle(v_x-\langle v_x\rangle)^2\rangle\right],
\end{eqnarray}
where $y$ is the lateral direction ($v_y=0$ in the two-dimensional simulations). The term on the right-hand side ($\langle(v_x-\langle v_x\rangle)^2\rangle$) represents the leading-order correlation term. Higher-order correlations and fluctuations of $u_x$ have been neglected. For simplicity, we assume that particle velocity fluctuations are isotropic, thus $\langle(v_x-\langle v_x\rangle)^2\rangle=\langle v_y^2\rangle=\langle v_z^2\rangle$ ($\langle v_y\rangle=\langle v_z\rangle=0$ for steady transport \cite{Paehtzetal15a}). Averaging over height approximately leads to
\begin{eqnarray}
 \frac{\overline{a^d_x}|_\mathrm{turb}}{\tilde g} = \frac{3}{4}C_{\rm fluc}C_d^\infty(U_x-V_x)^2,
\end{eqnarray}
where $C_{\rm fluc} = \sqrt{1+2V_z^2/(U_x-V_x)^2}+V_z^2/(U_x-V_x)^2$.
%Note that this drag correction applies to any constant drag coefficient and thus also to the coefficient `$2$' in Eq.~(\ref{Dragrare}).

\subsubsection*{Cohesion correction ($\mu_b(C)$ in Eq.~(1))}

We assume that the effect of cohesive forces on the dynamic friction coefficient $\mu_b$ can be modelled in a first approximation by the cohesion correction of the static threshold derived by Ref.~\cite{ClaudinAndreotti06} 
%
\begin{equation}
	\mu_b(C) = \mu_b^o\left[1+1.5(c_\mathrm{coh}C)^{5/3}\right], \label{mub2} 
\end{equation}
%
where $C$ is the cohesion number (defined in the text), $c_\mathrm{coh}$ a dimensionless constant characterizing the strength of cohesive forces, and $\mu_b^o$ the cohesionless friction coefficient. We obtain $c_\mathrm{coh}=0.19$ (similar to $c_\mathrm{coh}=0.11$ by Ref.~\cite{ClaudinAndreotti06}) from imposing that the smallest value of the threshold fluid shear velocity $u^\mathrm{ex}_t$ as a function of the particle diameter $d$ for conditions in Earth's atmosphere is predicted for $d=125\mu$m, consistent with Chepil's \cite{Chepil45} measurements, which exhibited a minimal value for $d=125\mu$m.

\subsubsection*{The law of the wall (function $f_2$ in Eq.~(2))}
The horizontal fluid velocity $u_x(z)$ at the threshold is calculated using the ``law of the wall'' proposed by Ref.~\cite{GuoJulien07},
\begin{eqnarray}
 \frac{u_x}{\sqrt{(s-1)gd}} & = & \sqrt{\Theta^{\mathrm{ex}}_t}f_2\left( (z-z_s) d^{-1} / Z_\nu \right), \nonumber \\
 f_2\left( z d^{-1} / Z_\nu \right) & = & 7\arctan\left( \frac{1}{7 Z_\nu}\frac{z+z_u}{d} \right)+\frac{7}{3}\arctan^3 \left( \frac{1}{7 Z_\nu}\frac{z+z_u}{d} \right) \nonumber \\ 
 && - 0.52\arctan^4\left( \frac{1}{7 Z_\nu}\frac{z+z_u}{d} \right) + \ln\left[1+\left(\frac{1}{B Z_\nu}\frac{z+z_u}{d}\right)^{(1/\kappa)}\right] \nonumber \\
 &&-\frac{1}{\kappa}\ln\left\{1+\frac{0.3}{Z_\nu} \left[1-\exp\left(-\frac{1}{26 Z_\nu}\right)\right]\right\}, \label{uxcomplex}
\end{eqnarray}
where $\kappa=0.4$, $B=\exp(16.873\kappa-\ln9)$, and $z_u\approx0.7d$ (from our direct transport simulations) takes into account that the fluid velocity profile extrapolated into the sediment bed vanishes a constant distance below $z_s$. This version has the advantage of providing a single equation for all regimes. Within the viscous sublayer of the turbulent boundary layer $u_x/\sqrt{(s-1)gd} \rightarrow \sqrt{\Theta^{\mathrm{ex}}_t}(z-z_u) d^{-1}/Z_\nu$, whereas in the log-layer $u_x/\sqrt{(s-1)gd} \rightarrow \kappa^{-1}\sqrt{\Theta^{\mathrm{ex}}_t}\ln((z-z_u)/z_o)$. The roughness length $z_o$ equals $Z_\nu d/9$ in the hydraulic smooth and $d/30$ in the hydraulic rough regime \cite{GuoJulien07}.

\subsubsection*{Derivation of Eq.~4}
From the vertical momentum balance \cite{Paehtzetal15a} $\d P_{zz}/\d z=-\rho\tilde g$, where $\rho$ is the particle mass density, and partial integration, we obtain
\begin{eqnarray}
Z = (\overline{z}-z_s) d^{-1} = \frac{\int_{z_s}^\infty\rho z\d z}{\int_{z_s}^\infty\rho d\d z} - z_s d^{-1} = \frac{\int_{z_s}^\infty\rho(P_{zz}/\rho)\d z}{P_{zz}d}=\frac{\overline{P_{zz}/\rho}}{\tilde g d}.
\end{eqnarray}
Using Eq.~(\ref{Particlestress}), we finally get
\begin{equation}
	Z =  Z_c + s V_z^2,
\end{equation}
where $V_z^2 = \overline{v_z^2}/[(s-1)gd]$ and $Z_c = \overline{P^c_{zz}/\rho} (\hat gd)^{-1}$ is proportional to $\mu_b^{-1}\Theta$ in our direct transport simulations (Extended Data Fig.~4).

\subsubsection*{Derivation of Eq.~5}
Using Eq.~(\ref{Particlestress}) and neglecting the contributions from contact forces $P^c_{ij}$ to the particle stress tensor, we can approximate the friction coefficient as \cite{Paehtzetal15a}
\begin{eqnarray}
 \mu = -\frac{P_{zx}}{P_{zz}}\approx-\frac{P^t_{zx}}{P^t_{zz}}=-\frac{\langle v_zv_x\rangle}{\langle v_z^2\rangle}. \label{muapprox}
\end{eqnarray}
Note that neglecting $P^c_{ij}$ is actually not justified for conditions with transport layer size of the order of $d$ (`bedload', such as subaqueous sediment transport) in our simulations, for which $P^c_{ij}$ constitutes the dominant contribution to the particle stress tensor. However, the error of the approximation in Eq.~(\ref{muapprox}) is nonetheless small (typically less than $30\%$) for these cases because the errors in the approximations of the numerator ($P_{zx}\approx\rho\langle v_zv_x\rangle$) and denominator ($P_{zz}\approx\rho\langle v_z^2\rangle$) of Eq.~(\ref{muapprox}) seem to largely compensate each other at locations above the bed surface for reasons we do not fully understand (not shown). Moreover, the influence of Eq.~(5) on the threshold calculation is very small for bedload conditions with $\mathrm{Re}_p\lessapprox50$, since the term $Z_c$ in Eq.~(4) dominates the term containing $V_z$ for bedload conditions, and since the contribution of $V_z$ in Eq.~(1) is only significant for bedload with $\mathrm{Re}_p\gtrapprox50$ (not shown). That is, neglecting interparticle contacts in the derivation of Eq.~(5) is justified in most cases in our simulations for the purpose of calculating the threshold.

Now, after averaging Eq.~(\ref{muapprox}) over height and assuming $\overline{v_z v_x}=\overline{|v_z|}(\overline{v_{x\downarrow}}-\overline{v_{x\uparrow}})\approx \alpha \overline{v_x}\overline{|v_z|}$, we obtain
\begin{eqnarray}
 \overline{\mu} \approx \alpha\frac{\overline{v_x}}{\overline{|v_z|}} = \alpha \frac{V_x}{\overline{|v_z|}/\sqrt{(s-1)gd}}. \label{beforeclosure}
\end{eqnarray}
Finally, we apply the closure $\overline{\mu}\approx\mu_b$ and $\overline{|v_z|}\approx\sqrt{\overline{v_z^2}}$, which yields
\begin{eqnarray}
 V_z=\alpha\mu_b^{-1}V_x. \label{afterclosure}
\end{eqnarray}
Although this closure significantly worsens the agreement with the simulation data (Extended Data Figures~5a versus 5b), it has the advantage of reducing the number of variables and thus model equations. Note that the simulations data indicate that $V_x$ in Eqs.~(\ref{beforeclosure}) and (\ref{afterclosure}) should be replaced by $\sqrt{V_x|_{\Theta=\Theta^{\mathrm{ex}}_t}V_x}$ for reasons we do not fully understand. However, this does not matter for the threshold calculation as both terms are the same near threshold conditions.

\subsubsection*{Characteristic travel height and length of fast particles}

From Eq.~(5), assuming ballistic trajectories, and using that the transport-layer-averaged particle velocity is a measure for the characteristic speed of fast transported particles \cite{PaehtzDuran16a}, the characteristic travel height $H$ and length $L$ of fast particles are given by $H/d=s V_z^2$ and $L/d=2 s V_x V_z = 2\alpha^{-1}\mu_b(C)H \gg H$, where $V_x$ and $V_z$ are defined in the text.

\section*{Supplementary Online Text}

\subsubsection*{1. Influence of particle hindrance}
The hindrance effect, which is neglected in our analytical model, increases the viscosity and thus the viscous drag force acting on particles in the presence of neighbouring particles by about factor $(1-\phi/\phi_{\mathrm{max}})^{-2.5\phi_{\mathrm{max}}}$ \cite{IshiiHibiki11}, where $\phi=\rho/\rho_p$ is the particle volume fraction and $\phi_{\mathrm{max}}\approx0.58$ its maximal value (random close packing). If the transport layer is sufficiently large, $\phi$ will be very small at typical transport heights, and this effect can thus be neglected. However, for viscous bedload (i.e., transport layer size of the order of $d$), the hindrance effect is likely significant since the transport layer is very thin ($Z\approxeq Z_c$), corresponding to comparably large values of $\phi$ within the transport layer. An average volume fraction of about $0.16$ would explain the amount of overestimation of $\Theta^{\mathrm{ex}}_t$ by our model in comparison to the measurements by Ref.~\cite{Ouremietal07}. This volume fraction is, indeed, a typical value within the transport layer of viscous bedload in our transport simulations, which is why we believe that the overestimation of the viscous measurements is mainly due to neglecting the hindrance effect.

\subsubsection*{2. The combined effect of the viscous sublayer and impact entrainment on the threshold}
It is instructive to analyse the analytical model predictions for the following ideal cases: (1) viscous saltation in non-rarefied air, (2) aerodynamically smooth turbulent saltation in non-rarefied air, for which the transport-layer-averaged fluid speed is predominantly controlled by the viscous sublayer and the average particle speed by the turbulent boundary layer, and (3) aerodynamically rough turbulent saltation in non-rarefied air, for which both speeds are controlled by the turbulent boundary layer. 

\subsubsection*{2.1 Viscous saltation in non-rarefied air}
For viscous saltation, we obtain $U_x-V_x=\mu_b\mathrm{Re}_p/18$ from Eq.~(1), $U_x=\Theta^{\mathrm{ex}}_t\mathrm{Re}_pZ$ from Eq.~(2), $V_x=\beta_1\sqrt{\beta_2}U_x$ from Eq.~(3), and $Z\approxeq\alpha^2\mu_b^{-2}sV_x^2$ (large $s$) from Eqs.~(4) and (5). Combined, this yields
\begin{eqnarray}
 \Theta^{\mathrm{ex}}_t\approxeq\frac{18\mu_b(1-\beta_1\sqrt{\beta_2})}{\alpha^2\beta_1^2\beta_2}(\mathrm{Im})^{-2}, \label{Visc}
\end{eqnarray}
where $\mathrm{Im}=\mathrm{Re}_p\sqrt{s+0.5}$ is the impact number \cite{PaehtzDuran16a}. Eq.~(\ref{Visc}) approximately implies $u^{\mathrm{ex}}_t\propto\mathrm{Re}_p^{-1}\sqrt{gd}$ (independent of $s$). Hence, the increase of the transport layer height $Z$ with $s$, which tends to increase the average fluid drag force acting on particles due to larger fluid velocities, exactly compensates the natural decrease of the average fluid drag force with $s$. The reason is the linear (=strong) increase of the fluid velocity with height within the viscous sublayer.

\subsubsection*{2.2 Aerodynamically smooth turbulent saltation in non-rarafied air}
However, the transport height of particles cannot increase indefinitely. At some point, the fastest moving particles will exceed the viscous sublayer and become influenced by the turbulent boundary layer, though fluid drag and the transport-layer-averaged fluid speed tend to be in the viscous regime. Then $V_x$ will become controlled by $\sqrt{\Theta^{\mathrm{ex}}_t}$ rather than $U_x$. This is encoded in the saturated value of Eq.~(3) (large $U_x/\sqrt{\Theta^{\mathrm{ex}}_t}$), which reads $V_x=\beta_1\sqrt{\Theta^{\mathrm{ex}}_t}$. In the saturated regime, $U_x/V_x$ will increase strongly as $U_x$ remains to be dominated by the viscous sublayer. We thus further approximate $U_x-V_x\approxeq U_x$. These two changes imply
\begin{eqnarray}
 \Theta^{\mathrm{ex}}_t\approxeq\sqrt{\frac{\mu_b^3}{18\alpha^2\beta_1^2}}\sqrt{s}^{-1},
\end{eqnarray}
which approximately implies $u^{\mathrm{ex}}_t\propto s^{1/4}\sqrt{gd}$ (independent of $\mathrm{Re}_p$). Aerodynamically smooth turbulent saltation in non-rarafied air is the most relevant regime on Earth and Mars, and also relevant on Triton and Pluto.

\subsubsection*{2.3 Aerodynamically rough turbulent saltation in non-rarafied air}
For aerodynamically rough turbulent saltation, the viscous sublayer becomes negligible, and fluid drag and the transport-layer-averaged fluid speed tend to be in the turbulent regime. We thus obtain $U_x-V_x\approxeq U_x=\sqrt{4\mu_b/(3C_d^\infty)}$ from Eq.~(1), $U_x=\kappa^{-1}\sqrt{\Theta^{\mathrm{ex}}_t}\ln(Zd/z_o)$ from Eq.~(2), $V_x=\beta_1\sqrt{\Theta^{\mathrm{ex}}_t}$ from Eq.~(3), and $Z\approxeq\alpha^2\mu_b^{-2}sV_x^2$ (large $s$) from Eqs.~(4) and (5). Combined, this yields
\begin{eqnarray}
 \sqrt{\Theta^{\mathrm{ex}}_t}=\frac{\kappa\sqrt{4\mu_b}}{\sqrt{3C_d^\infty}\ln\left(\frac{\alpha^2\Theta^{\mathrm{ex}}_tsd}{\mu_b^2z_o}\right)}, \label{Turb}
\end{eqnarray}
which approximately implies $u^{\mathrm{ex}}_t\propto (\ln s)^{-1}s^{1/2}\sqrt{gd}$. One can see that now, when the viscous sublayer does not play an important role, $u^{\mathrm{ex}}_t$ increases much more strongly with $s$ than in the previous two cases, which outlines the immense importance of the viscous sublayer for low threshold values when $s$ is large (e.g., on Triton, Pluto, comet 67P).

The predicted transport regimes (Eqs.~(\ref{Visc}-\ref{Turb})) are outlined in Extended Data Figure~7, which shows predicted values of $\Theta^{\mathrm{ex}}_t$ (solid lines) for non-rarefied ($\mathrm{Kn}=0$), non-cohesive ($C=0$) conditions as a function of the impact number ($\mathrm{Im}$) for different values of $s$. It also shows the predictions of the model of Ref.~\cite{ClaudinAndreotti06} (dashed lines), which are quantitatively and qualitatively different because the effects of near-bed collisions were modelled differently.


\section*{Movies}
\subsection{Movie S1 - Particle-bed impact entrainment in subaqueous sediment transport}
Time evolution of the simulated particle-fluid system for $s=2.65$, $\mathrm{Re}_p=20$, $\sqrt{\Theta}=0.33$ (near threshold conditions), and nearly maximally damped interparticle collisions (restitution coefficient $e=0.01$). The flow velocity is shown as a background colour with warm colours corresponding to high velocities and cold colours to small velocities. The horizontal and vertical axes are measured in mean particle diameters. Only $1/4$ of the simulated horizontal domain is shown, which is why particles occasionally enter the system from the left. It can be seen that impacting particles occasionally drag bed particles out of their traps.

\subsection{Movie S2 - Particle-bed impact entrainment in aeolian sediment transport}
Time evolution of the simulated particle-fluid system for $s=2000$, $\mathrm{Re}_p=20$, $\sqrt{\Theta}=0.08$ (near threshold conditions), and weakly damped interparticle collisions (restitution coefficient $e=0.9$). The flow velocity is shown as a background colour with warm colours corresponding to high velocities and cold colours to small velocities. The horizontal and vertical axes are measured in mean particle diameters. Only $1/4$ of the simulated horizontal domain is shown, which is why particles occasionally enter the system from the left. This is an example for sediment transport sustained predominantly through particle-bed-impact entrainment. It can be seen that impacting particles occasionally eject bed particles.

\end{document}

1602.07079.extracted_bib
{Hansenetal90 \bibinfo{author{Hansen, C.~J., \bibinfo{author{McEwen, A.~S., \bibinfo{author{Ingersoll, A.~P. \& \bibinfo{author{Terrile, R.~J. \newblock \bibinfo{title{Surface and airborne evidence for plumes and winds on triton. \newblock \bibinfo{journal{Science \bibinfo{volume{250, \bibinfo{pages{421--424 (\bibinfo{year{1990).
1602.07079.extracted_bib
{SaganChyba90 \bibinfo{author{Sagan, C. \& \bibinfo{author{Chyba, C. \newblock \bibinfo{title{Triton's streaks as windblown dust. \newblock \bibinfo{journal{Nature \bibinfo{volume{346, \bibinfo{pages{546--548 (\bibinfo{year{1990).
1602.07079.extracted_bib
{Sternetal15 \bibinfo{author{Stern, S.~A. et~al. \newblock \bibinfo{title{The pluto system: Initial results from its exploration by new horizons. \newblock \bibinfo{journal{Science \bibinfo{volume{350, \bibinfo{pages{aad1815 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Thomasetal15a \bibinfo{author{Thomas, N. et~al. \newblock \bibinfo{title{The morphological diversity of comet 67p/churyumov-gerasimenko. \newblock \bibinfo{journal{Science \bibinfo{volume{347, \bibinfo{pages{aaa0440 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Thomasetal15b \bibinfo{author{Thomas, N. et~al. \newblock \bibinfo{title{Redistribution of particles across the nucleus of comet 67p/churyumov-gerasimenko. \newblock \bibinfo{journal{Astronomy \& Astrophysics \bibinfo{volume{583, \bibinfo{pages{A17 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Mooreetal15 \bibinfo{author{Moore, J.~M. et~al. \newblock \bibinfo{title{Geology before pluto: Pre-encounter considerations. \newblock \bibinfo{journal{Icarus \bibinfo{volume{246, \bibinfo{pages{65--81 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Gladstoneetal16 \bibinfo{author{Gladstone, G.~R. et~al. \newblock \bibinfo{title{The atmosphere of pluto as observed by new horizons. \newblock \bibinfo{journal{Science \bibinfo{volume{351, \bibinfo{pages{aad8866 (\bibinfo{year{2016).
1602.07079.extracted_bib
{Kok10a \bibinfo{author{Kok, J.~F. \newblock \bibinfo{title{Difference in the wind speeds required for initiation versus continuation of sand transport on mars: Implications for dunes and dust storms. \newblock \bibinfo{journal{Physical Review Letters \bibinfo{volume{104, \bibinfo{pages{074502 (\bibinfo{year{2010).
1602.07079.extracted_bib
{ClaudinAndreotti06 \bibinfo{author{Claudin, P. \& \bibinfo{author{Andreotti, B. \newblock \bibinfo{title{A scaling law for aeolian dunes on mars, venus, earth, and for subaqueous ripples. \newblock \bibinfo{journal{Earth and Planetary Science Letters \bibinfo{volume{252, \bibinfo{pages{30--44 (\bibinfo{year{2006).
1602.07079.extracted_bib
{Kok10b \bibinfo{author{Kok, J.~F. \newblock \bibinfo{title{An improved parameterization of wind blown sand flux on mars that includes the effect of hysteresis. \newblock \bibinfo{journal{Geophysical Research Letters \bibinfo{volume{37, \bibinfo{pages{L12202 (\bibinfo{year{2010).
1602.07079.extracted_bib
{Paehtzetal12 \bibinfo{author{P\"ahtz, T., \bibinfo{author{Kok, J.~F. \& \bibinfo{author{Herrmann, H.~J. \newblock \bibinfo{title{The apparent roughness of a sand surface blown by wind from an analytical model of saltation. \newblock \bibinfo{journal{New Journal of Physics \bibinfo{volume{14, \bibinfo{pages{043035 (\bibinfo{year{2012).
1602.07079.extracted_bib
{Duranetal12 \bibinfo{author{Dur\'an, O., \bibinfo{author{Andreotti, B. \& \bibinfo{author{Claudin, P. \newblock \bibinfo{title{Numerical simulation of turbulent sediment transport, from bed load to saltation. \newblock \bibinfo{journal{Physics of Fluids \bibinfo{volume{24, \bibinfo{pages{103306 (\bibinfo{year{2012).
1602.07079.extracted_bib
{PaehtzDuran16a \bibinfo{author{P\"ahtz, T. \& \bibinfo{author{Dur\'an, O. \newblock \bibinfo{title{Fluid forces or impacts, what governs the entrainment of soil particles in sediment transport mediated by a newtonian fluid? \newblock \bibinfo{journal{under review in Physical Review Letters (\bibinfo{year{2016).
1602.07079.extracted_bib
{Clarketal15 \bibinfo{author{Clark, A.~H., \bibinfo{author{Shattuck, M.~D., \bibinfo{author{Ouellette, N.~T. \& \bibinfo{author{O'Hern, C.~S. \newblock \bibinfo{title{Onset and cessation of motion in hydrodynamically sheared granular beds. \newblock \bibinfo{journal{Physical Review E \bibinfo{volume{92, \bibinfo{pages{042202 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Maurinetal15 \bibinfo{author{Maurin, R., \bibinfo{author{Chauchat, J., \bibinfo{author{Chareyre, B. \& \bibinfo{author{Frey, P. \newblock \bibinfo{title{A minimal coupled fluid-discrete element model for bedload transport. \newblock \bibinfo{journal{Physics of Fluids \bibinfo{volume{27, \bibinfo{pages{113302 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Koketal12 \bibinfo{author{Kok, J.~F., \bibinfo{author{Parteli, E. J.~R., \bibinfo{author{Michaels, T.~I. \& \bibinfo{author{Karam, D.~B. \newblock \bibinfo{title{The physics of wind-blown sand and dust. \newblock \bibinfo{journal{Reports on Progress in Physics \bibinfo{volume{75, \bibinfo{pages{106901 (\bibinfo{year{2012).
1602.07079.extracted_bib
{Bagnold56 \bibinfo{author{Bagnold, R.~A. \newblock \bibinfo{title{The flow of cohesionless grains in fluid. \newblock \bibinfo{journal{Philosophical Transactions of the Royal Society London A \bibinfo{volume{249, \bibinfo{pages{235--297 (\bibinfo{year{1956).
1602.07079.extracted_bib
{Camenen07 \bibinfo{author{Camenen, B. \newblock \bibinfo{title{Simple and general formula for the settling velocity of particles. \newblock \bibinfo{journal{Journal of Hydraulic Engineering \bibinfo{volume{133, \bibinfo{pages{229--233 (\bibinfo{year{2007).
1602.07079.extracted_bib
{Croweetal12 \bibinfo{author{Crowe, C.~T., \bibinfo{author{Schwarzkopf, J.~D., \bibinfo{author{Sommerfeld, M. \& \bibinfo{author{Tsuji, Y. \newblock \bibinfo{title{Multiphase Flows with Droplets and Particles (\bibinfo{publisher{Taylor \& Francis Group, Boca Raton, \bibinfo{year{2012).
1602.07079.extracted_bib
{Buffingtonmontgomery97 \bibinfo{author{Buffington, J.~M. \& \bibinfo{author{Montgomery, D.~R. \newblock \bibinfo{title{A systematic analysis of eight decades of incipient motion studies, with special reference to gravel-bedded rivers. \newblock \bibinfo{journal{Water Resources Research \bibinfo{volume{33, \bibinfo{pages{1993--2029 (\bibinfo{year{1997).
1602.07079.extracted_bib
{Ouremietal07 \bibinfo{author{Ouriemi, M., \bibinfo{author{Aussillous, P., \bibinfo{author{Medale, M., \bibinfo{author{Peysson, Y. \& \bibinfo{author{Guazzelli, E. \newblock \bibinfo{title{Determination of the critical shields number for particle erosion in laminar flow. \newblock \bibinfo{journal{Physics of Fluids \bibinfo{volume{19, \bibinfo{pages{061706 (\bibinfo{year{2007).
1602.07079.extracted_bib
{Bagnold37 \bibinfo{author{Bagnold, R.~A. \newblock \bibinfo{title{The transport of sand by wind. \newblock \bibinfo{journal{The Geographical Journal \bibinfo{volume{89, \bibinfo{pages{409--438 (\bibinfo{year{1937).
1602.07079.extracted_bib
{Chepil45 \bibinfo{author{Chepil, W.~S. \newblock \bibinfo{title{Dynamics of wind erosion: Ii. initiation of soil movement. \newblock \bibinfo{journal{Soil Science \bibinfo{volume{60, \bibinfo{pages{397--411 (\bibinfo{year{1945).
1602.07079.extracted_bib
{Creysselsetal09 \bibinfo{author{Creyssels, M. et~al. \newblock \bibinfo{title{Saltating particles in a turbulent boundary layer: experiment and theory. \newblock \bibinfo{journal{Journal of Fluid Mechanics \bibinfo{volume{625, \bibinfo{pages{47--74 (\bibinfo{year{2009).
1602.07079.extracted_bib
{Hoetal11 \bibinfo{author{Ho, T.~D., \bibinfo{author{Valance, A., \bibinfo{author{Dupont, P. \& \bibinfo{author{{Ould El Moctar, A. \newblock \bibinfo{title{Scaling laws in aeolian sand transport. \newblock \bibinfo{journal{Physical Review Letters \bibinfo{volume{106, \bibinfo{pages{094501 (\bibinfo{year{2011).
1602.07079.extracted_bib
{Sternetal16 \bibinfo{author{Stern, S.~A. et~al. \newblock \bibinfo{title{Evidence that pluto's long-term polar axis shifts drive sharp changes in the planet’s atmospheric pressure over time. \newblock \bibinfo{journal{under review in Icarus (\bibinfo{year{2016). \newblock \bibinfo{note{Video (`Press Briefing - New Horizons: Peering into Pluto’s Past'): http://livestream.com/viewnow/LPSC2016.
1602.07079.extracted_bib
{Trilling16 \bibinfo{author{Trilling, D.~E. \newblock \bibinfo{title{The surface age of sputnik planum, pluto, must be less than 10 million years. \newblock \bibinfo{journal{PLOS ONE \bibinfo{volume{11, \bibinfo{pages{e0147386 (\bibinfo{year{2016).
1602.07079.extracted_bib
{Hand15 \bibinfo{author{Hand, E. \newblock \bibinfo{title{Late harvest from pluto reveals a complex world. \newblock \bibinfo{journal{Science \bibinfo{volume{350, \bibinfo{pages{260--261 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Burretal15 \bibinfo{author{Burr, D.~M. et~al. \newblock \bibinfo{title{Higher-than-predicted saltation threshold wind speeds on titan. \newblock \bibinfo{journal{Nature \bibinfo{volume{517, \bibinfo{pages{60--63 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Ayoubetal14 \bibinfo{author{Ayoub, F. et~al. \newblock \bibinfo{title{Threshold for sand mobility on mars calibrated from seasonal variations of sand flux. \newblock \bibinfo{journal{Nature Communication \bibinfo{volume{5, \bibinfo{pages{5096 (\bibinfo{year{2014).
1602.07079.extracted_bib
{Tokano10 \bibinfo{author{Tokano, T. \newblock \bibinfo{title{Relevance of fast westerlies at equinox for the eastward elongation of titan's dunes. \newblock \bibinfo{journal{Aeolian Research \bibinfo{volume{2, \bibinfo{pages{113--127 (\bibinfo{year{2010).
1602.07079.extracted_bib
{Fentonetal05 \bibinfo{author{Fenton, L.~K., \bibinfo{author{Toigo, A.~D. \& \bibinfo{author{Richardson, M.~I. \newblock \bibinfo{title{Aeolian processes in proctor crater on mars: Mesoscale modeling of dune-forming winds. \newblock \bibinfo{journal{Journal of Geophysical Research \bibinfo{volume{110, \bibinfo{pages{E06005 (\bibinfo{year{2005).
1602.07079.extracted_bib
{Bagnold38 \bibinfo{author{Bagnold, R.~A. \newblock \bibinfo{title{The measurement of sand storms. \newblock \bibinfo{journal{Proceedings of the Royal Society London Series A \bibinfo{volume{167, \bibinfo{pages{282--291 (\bibinfo{year{1938).
1602.07079.extracted_bib
{GreeleyArvidson90 \bibinfo{author{Greeley, R. \& \bibinfo{author{Arvidson, R.~E. \newblock \bibinfo{title{Aeolian processes on venus. \newblock \bibinfo{journal{Earth, Moon, and Planets \bibinfo{volume{50, \bibinfo{pages{127--157 (\bibinfo{year{1990).
1602.07079.extracted_bib
{Greenstreetetal15 \bibinfo{author{Greenstreet, S., \bibinfo{author{Gladman, B. \& \bibinfo{author{McKinnon, W.~B. \newblock \bibinfo{title{Impact and cratering rates onto pluto. \newblock \bibinfo{journal{Icarus \bibinfo{volume{258, \bibinfo{pages{267--288 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Zahnleetal03 \bibinfo{author{Zahnle, K., \bibinfo{author{Schenk, P., \bibinfo{author{Levison, H. \& \bibinfo{author{Donesc, L. \newblock \bibinfo{title{Cratering rates in the outer solar system. \newblock \bibinfo{journal{Icarus \bibinfo{volume{163, \bibinfo{pages{263--289 (\bibinfo{year{2003).
1602.07079.extracted_bib
{FlynnMcKay90 \bibinfo{author{Flynn, G.~J. \& \bibinfo{author{McKay, D.~S. \newblock \bibinfo{title{An assessment of the meteoritic contribution to the martian soil. \newblock \bibinfo{journal{Journal of Geophysical Research \bibinfo{volume{95, \bibinfo{pages{14497--14509 (\bibinfo{year{1990).
1602.07079.extracted_bib
{McKinnonetal08 \bibinfo{author{McKinnon, W.~B., \bibinfo{author{Prialnik, D., \bibinfo{author{Stern, S.~A. \& \bibinfo{author{Coradini, A. \newblock \bibinfo{title{Structure and evolution of kuiper belt objects and dwarf planets. \newblock In \bibinfo{editor{Barucci, M.~A., \bibinfo{editor{Boehnhardt, H., \bibinfo{editor{Cruikshank, D.~P. \& \bibinfo{editor{Morbidelli, A. (eds.) \bibinfo{booktitle{The Solar System Beyond Neptune, \bibinfo{pages{213--241 (\bibinfo{publisher{University of Arizona Press, Tucson, \bibinfo{year{2008).
1602.07079.extracted_bib
{Paehtzetal15a \bibinfo{author{P\"ahtz, T., \bibinfo{author{Dur\'an, O., \bibinfo{author{Ho, T.-D., \bibinfo{author{Valance, A. \& \bibinfo{author{Kok, J.~F. \newblock \bibinfo{title{The fluctuation energy balance in non-suspended fluid-mediated particle transport. \newblock \bibinfo{journal{Physics of Fluids \bibinfo{volume{27, \bibinfo{pages{013303 (\bibinfo{year{2015).
1602.07079.extracted_bib
{Carneiroetal11 \bibinfo{author{Carneiro, M.~V., \bibinfo{author{P\"ahtz, T. \& \bibinfo{author{Herrmann, H.~J. \newblock \bibinfo{title{Jump at the onset of saltation. \newblock \bibinfo{journal{Physical Review Letters \bibinfo{volume{107, \bibinfo{pages{098001 (\bibinfo{year{2011).
1602.07079.extracted_bib
{GuoJulien07 \bibinfo{author{Guo, J. \& \bibinfo{author{Julien, P.~Y. \newblock \bibinfo{title{Buffer law and transitional roughness effect in turbulent open-channel flows. \newblock In \bibinfo{booktitle{The Fifth International Symposium on Environmental Hydraulics (ISEH V), Tempe, Arizona, vol. \bibinfo{volume{Paper 5, \bibinfo{pages{1--6 (\bibinfo{publisher{University of Nebraska - Lincoln, \bibinfo{year{2007).
1602.07079.extracted_bib
{IshiiHibiki11 \bibinfo{author{Ishii, M. \& \bibinfo{author{Hibiki, T. \newblock \bibinfo{title{of Two-Phase Flow Thermo-Fluid Dynamics, 323 pp. (\bibinfo{publisher{Springer, New York, \bibinfo{year{2011).