# Hydrodynamic roughness induced by a multiscale topography

Pan Jia1<EMAIL_ADDRESS>Bruno Andreotti2  Philippe Claudin3  1School of
Science, Harbin Institute of Technology, 518055 Shenzhen, P. R. China
2Laboratoire de Physique de l’ENS, UMR 8550 Ecole Normale Supérieure – CNRS –
Université PSL – Université Paris Cité – Sorbonne Université, 24 rue Lhomond,
75005 Paris, France 3Physique et Mécanique des Milieux Hétérogènes, UMR 7636
CNRS – ESPCI Paris – Université PSL – Université Paris Cité – Sorbonne
Université, 10 rue Vauquelin, 75005 Paris, France

###### Abstract

Turbulent flows above a solid surface are characterised by a hydrodynamic
roughness that represents, for the far velocity field, the typical length
scale at which momentum mixing occurs close to the surface. Here, we are
theoretically interested in the hydrodynamic roughness induced by a two-
dimensional modulated surface, the elevation profile of which is decomposed in
Fourier modes. We describe the flow for a sinusoidal mode of given wavelength
and amplitude with RANS equations closed by means of a mixing-length approach
that takes into account a possible surface geometrical roughness as well as
the presence of a viscous sublayer. It also incorporates spatial transient
effects at the laminar-turbulent transition. Performing a weekly non-linear
expansion in the bedform aspect ratio, we predict the effective hydrodynamic
roughness when the surface wavelength is varied and we show that it presents a
non-monotonic behaviour at the laminar-turbulent transition when the surface
is hydrodynamically smooth. Further, with a self-consistent looped
calculation, we are able to recover the smooth-rough transition of a flat
surface, for which the hydrodynamic roughness changes from a regime where it
is dominated by the viscous length to another one where it scales with the
surface corrugation. We finally apply the results to natural patterns resulted
from hydrodynamic instabilities such as those associated with dissolution or
sediment transport. We discuss in particular the aspect ratio selection of
dissolution bedforms and roughness hierarchy in superimposed ripples and
dunes.

###### keywords:

##  1 Introduction

The length scale over which mixing occurs in the superficial layer of a
turbulent boundary layer, close to a solid surface, is a key quantity in
meteorology and geomorphology. This length, later defined as the ‘hydrodynamic
roughness’, is for example necessary to account for the effect of relief and
vegetation on winds (Gillies et al., 2007; King et al., 2006; Lancaster &
Baas, 1998; Marticorena & Bergametti, 1995; Raupach, 1992; Raupach et al.,
1993; Wolfe & Nickling, 1993; Finnigan, 1988; Bradley, 1980). Remote sensing
methods including satellite imagery (Jasinski & Crago, 1999), airborne LIDAR
(Paul-Limoges et al., 2013), tethersonde and eddy correlation systems (Tsai et
al., 2010; Tsuang et al., 2003), GPS radiosondes (Han et al., 2015), ERS
scattermeter (Prigent et al., 2005), and terrestrial laser scanning (Nield et
al., 2013), have been employed to evaluate the hydrodynamic roughness and its
relation to geometrical asperities of the ground. In glaciology, these
measurements are also used to quantify the turbulent heat exchange between the
surface of glaciers and the atmosphere, and their contribution to the ice melt
(Smeets & Van den Broeke, 2008). Similarly, the determination of the
hydrodynamic roughness of the sea surface is important in oceanography to
quantify the exchange of momentum at the air-sea interface (Maat et al., 1991)
and the nucleation of waves (Perrard et al., 2019; Aulnette et al., 2019). As
a final example, this length scale directly controls the vertical thermal
exchange with the atmosphere in cities (Kent et al., 2017; Crago et al.,
2012).

Some of these applied problems belong to the generic case of a turbulent
boundary layer above a flat rough surface. Since the earliest investigations
in the 1930s, including those of Nikuradse (1933) and Schlichting (1937), it
has been found that the velocity profile in such flows is logarithmic with
respect to the distance $z$ to the surface. Following dimensional analysis
(Prandtl, 1925; Schlichting & Gersten, 2000; Pope, 2000), the rate of momentum
mixing in this turbulent layer is determined by the mean velocity gradient,
over a length scale, called the mixing length, $\ell$, proportional to the
geometrical distance $z$. For a steady and homogenous flow over a flat
substrate in neutral conditions, the velocity profile away from the boundary
takes the classical form:

$u_{x}=\frac{u_{*}}{\kappa}\ln\left(\frac{z}{z_{0}}\right)\,,$ (1)

where $u_{x}$ is the time-averaged velocity along the flow direction $x$.
$u_{*}$ is the friction velocity, related by $u_{*}=\sqrt{\tau_{xz}/\rho}$ to
the shear stress $\tau_{xz}$ and the fluid density $\rho$. $\tau_{xz}$
represents the flux of momentum across surfaces parallel to the solid
substrate and is therefore conserved in the $z$ direction. $\kappa\approx 0.4$
is the phenomenological von Kármán constant. The length $z_{0}$ is the
altitude at which the logarithmic velocity profile (1) seems to vanish above
the ground, when extrapolated, hence appearing as a _reference level_ below
which the velocity becomes very small. It is called the hydraulic roughness in
hydraulic engineering (Van Rijn, 1984; Flack & Schultz, 2010), and the
aerodynamic roughness in the case of an air flow, and here generically named
hydrodynamic roughness. $z_{0}$ can be interpreted as a mixing length
governing the turbulent fluctuations close to the solid surface. The physical
origin of the hydrodynamic roughness depends on the processes at play close to
the surface. It has for instance been experimentally studied in wind tunnels
with discrete roughness elements placed on a flat substrate (Taylor et al.,
1976; Schmid & Bünzli, 1995; Hobson et al., 1999; Sadique et al., 2017). It is
modified when the solid is covered by slender flexible structures such as
plants (De Langre, 2008). It can be produced by turbulence fluctuations
themselves, by roughening of the surface of the ocean. Another example is the
case of sediment transport, which can dominate the interaction between the
granular bed surface and the flow, and typically increases the hydrodynamic
roughness (Owen, 1964; Raupach, 1991; Gillette et al., 1998; Sherman &
Farrell, 2008; Durán et al., 2011).

Figure 1:  (a) Schematic of superimposed bedforms: grains of size $d$ form
sand ripples of wavelength $\lambda$ and amplitude $\zeta$, which themselves
cover the surface of a dune; the hydrodynamic roughness at the scale of the
large bedform is induced by the small ones. (b) Martian dune of Nili Patera
dune field (from HIRISE picture ESP 062069 1890) with superimposed meter-scale
ripples. Dune size $\simeq 300$ m. (c) Photograph of the surface of a Martian
dune in Gale Crater by NASA’s Curiosity rover (Image Credit: NASA/JPL-
Caltech/MSSS PIA20755). Small decimetric ripples superimposed on large meter-
scale bedforms.

A situation of interest is that of a turbulent flow over an elevation profile
that presents several nested length scales. This situation of a gentle
topography with superimposed patterns is typically relevant of sedimentary
bedforms (Van Rijn, 1982, 1984; Venditti et al., 2005; Elbelrhiti et al.,
2005; Narteau et al., 2009; Nield et al., 2013; Durán Vinent et al., 2019). A
generic example is grains that form sand ripples, themselves covering the
surface of dunes. As shown in figure 1, sedimentary bedforms can present up to
four interlocking scales: grains, small impact ripples, large hydrodynamic
ripples and dunes. How to account for the fact that the grains control the
hydrodynamic roughness for the flow over the sand ripples, and the ripples
that over the dunes? How to handle such a multi-scale roughness hierarchy?
Here, we consider a two-scale problem, and investigate the hydrodynamic
roughness induced by a large scale topography whose surface presents a
roughness at a much smaller scale. What is the effect of the _inner_ surface
roughness on the _outer_ topography-induced roughness, far from the surface?
What are the separate effects of amplitude and horizontal length scale of the
large scale elevation profile on the outer hydrodynamic roughness? Is there a
signature of the laminar-turbulent transition on this hydrodynamic roughness?

The article is organised as follows. In section 2, we perform the weakly
nonlinear analysis of the turbulent flow above a modulated surface using two
expansion techniques. We first discuss the smooth-rough transition, and
present the turbulent closure used to relate the Reynolds stresses to the
velocity gradient inspired from models developed for erosion, dissolution and
sublimation pattern formation (Charru et al., 2013; Claudin et al., 2017;
Durán Vinent et al., 2019). This closure involves a mixing length which allows
us not only to account for surface geometrical roughness and viscous sublayer,
but also to incorporate transient effects at the laminar-turbulent transition.
The most technical mathematical considerations are gathered in several
appendices. The results are presented in section 3 where we study the outer
hydrodynamic roughness induced by a topography, with an emphasis on the effect
of the inner, surface roughness. We show that a self-consistent, looped
calculation, allows us to recover the hydrodynamic roughness induced by a
rough surface of equivalent grain size $d$. We then come back to the questions
mentioned in the previous paragraph, compare some of our findings with salient
results from the literature, and address the laminar-turbulent transition from
the point of view of bed roughness. Open problems and perspectives conclude
the manuscript.

Figure 2: (a) Hydrodynamic roughness $z_{0}$ measurements, including Flack &
Schultz (2010) (green squares) and those by Nikuradse (1933) (violet
triangles), as a function of the dimensionless equivalent grain size
$\mathcal{R}_{d}=du_{*}/\nu$. Green background $\mathcal{R}_{d}\lesssim 10$:
smooth regime for which $z_{0}\simeq\nu/(7u_{*})$. Gray background
$\mathcal{R}_{d}\gtrsim 10^{2}$: rough regime for which $z_{0}\simeq d/30$.
Lines: theoretical prediction of an effective roughness $z_{e}$ induced by a
sinusoidal surface of wavelength $\lambda=d$ and amplitude $\zeta=0.88d/2\pi$
(see section 3.3). Dash-dotted orange line: perturbative calculation (Eq. 16).
Solid red line: self-consistent calculation. (b) Dimensionless velocity
profiles $u_{x}/u_{*}$ from data over a smooth flat surface Johnson & Johnston
(1989); Wei & Willmarth (1989) (inverted triangles) and over various rough
flat surfaces (Schultz & Flack, 2009; Flack & Schultz, 2010) (squares). Lines:
function $\mathcal{U}$ (Eq. 10) in the smooth (black curve,
$\mathcal{R}_{d}=0$) and rough (blue curve, $\mathcal{R}_{d}\to\infty$)
limits. The distance to the surface $z$ is rescaled by $z_{0}$, deduced by an
extrapolation to a vanishing value of the profile semi-log fit far from the
surface (dashed orange line).

##  2 Hydrodynamic model

###  2.1 Smooth rough transition

We consider here a turbulent boundary layer above a solid surface whose
elevation profile presents two well-separated length scales: a large-scale
topography whose surface presents a roughness at a smaller scale. This
situation is typically that of sand grains at the surface of sand ripples, or
that of sand ripples at the surface of sand dunes. A simplifying description
is to consider the large scale topography as smooth, described by a vertical
elevation profile $z=Z(x)$, invariant along the spanwise $y$-direction, and to
take into account the effect of the roughness of the solid surface in the
hydrodynamic model (figure 1a).

The situation where the topography is flat provides important constrains to
calibrate such a model (figure 2). When the surface roughness is much smaller
than the thickness of the viscous sublayer $\sim\nu/u_{*}$, where $\nu$ is the
kinematic viscosity of the fluid, the hydrodynamic regime is said _smooth_ :
the logarithmic velocity profile is regularised by a laminar flow close to the
surface (Schlichting & Gersten, 2000; Pope, 2000). In this case, the
hydrodynamic roughness is consequently found to be controlled by the viscous
length $\nu/u_{*}$, independent of the amplitude of the solid surface
asperities and is around $z_{0}\sim\nu/(7u_{*})$ (Raupach et al., 1991b;
Schlichting & Gersten, 2000). In contrast, when the viscous sublayer is thin
enough, $z_{0}$ is rather controlled by the geometrical roughness of the solid
surface and the hydrodynamic regime is said _rough_. In this context, the
reference situation is a granular bed composed of natural sand grains of size
$d$, which leads to a hydrodynamic roughness $z_{0}=rd$ in the rough regime.
For a static sand bed of grain size $d$, the hydrodynamic roughness measured
by extrapolation of the velocity log-profile (1) in the rough regime gives
typical values between $r\simeq 0.03-0.1$ (Bagnold, 1941; Kamphuis, 1974;
Schlichting & Gersten, 2000; Andreotti, 2004). For an arbitrary rough
substrate it has been proposed to defined an equivalent sand grain such that
its hydrodynamic roughness is the same as that of a sand bed composed of
grains of that size (Flack & Schultz, 2010; Chung et al., 2021; Kadivar et
al., 2021). This approach assumes that the hydrodynamic roughness is
essentially proportional to the average roughness element height (Fang & Sill,
1992). Here, we will show, however, that the hydrodynamic roughness presents
more subtle dependences on other parameters such as the element spacing
density and shape (Wiberg & Nelson, 1992; Xue et al., 2002; Dong et al., 2002;
Jiménez, 2004; Cheng et al., 2007; Brown et al., 2008).

###  2.2 Governing equations and turbulent closure

In a spirit similar to Taylor et al. (1989), we consider that the inner
surface roughness is characterised by an equivalent grain size $d$. The
smooth-rough transition is controlled by the grain size based Reynolds number
$\mathcal{R}_{d}=du_{*}/\nu$. The velocity $\vec{u}$ and pressure $p$ fields
are described by means of Reynolds averaged Navier-Stokes equations, which
states

$\displaystyle\partial_{i}u_{i}$ $\displaystyle=$ $\displaystyle 0\,,$ (2)
$\displaystyle\rho u_{i}\partial_{i}u_{j}$ $\displaystyle=$
$\displaystyle-\partial_{j}p+\partial_{i}\tau_{ij}\,,$ (3)

where $\rho$ is the fluid density and $\tau_{ij}$ is the Reynolds stress
tensor. As a classical mixing theory for boundary layer flows, we take here a
first order closure (Pope, 2000) to relate the stress to the strain rate
$\dot{\gamma}_{ij}=\partial_{i}u_{j}+\partial_{j}u_{i}$, which adds
microscopic molecular and turbulent viscosities:

$\tau_{ij}=\rho(\nu+\ell^{2}|\dot{\gamma}|)\dot{\gamma}_{ij}-\frac{1}{3}\rho\chi^{2}\ell^{2}|\dot{\gamma}|^{2}\delta_{ij}\,,$
(4)

where
$|\dot{\gamma}|=\left(\tfrac{1}{2}\dot{\gamma}_{ij}\dot{\gamma}_{ij}\right)^{1/2}$
is the mixing frequency proportional to the modulus of the strain rate. $\chi$
is a phenomenological constant, which is in the range $2$–$3$ (Pope, 2000). It
has no importance here as only the normal stress difference
$\tau_{xx}-\tau_{zz}$ matters in the calculations (Fourrière et al., 2010;
Claudin et al., 2017). $\ell$ is the mixing length, which is the key input of
the model. The simplest expression of the mixing length presenting a smooth-
rough transition is $\ell=\kappa(z+rd)$, where the grain-induced roughness
$rd$ is added to the geometrical distance to the solid, $z$. However, this
formula overestimates turbulent mixing in the viscous sublayer, in the smooth
hydrodynamical regime. A classical empirical approach is to introduce a factor
to the geometrical length contribution (Van Driest, 1956) which exponentially
kills turbulent mixing when the local Reynolds number $u_{*}z/\nu$ is smaller
than a transitional value $\mathcal{R}^{0}_{t}\simeq 25$ (Pope, 2000). We
adopt here an expression introduced and tuned in the context of the linear
stability analysis for sedimentary or dissolution bedforms (Richards, 1980;
Ayotte et al., 1994; Colombini, 2004; Fourrière et al., 2010; Charru et al.,
2013; Claudin et al., 2017)

$\ell=\kappa(z+rd-Z)\left[1-\exp\left(-\frac{(\tau_{xz}/\rho)^{1/2}(z+sd-Z)}{\nu\mathcal{R}_{t}}\right)\right].$
(5)

$\kappa$ is the von Kármán constant as in Eq. 1. The parameter $s$ controls
the induction of turbulent fluctuations in the viscous layer upon increasing
the surface roughness. Both $r\simeq 1/30$ and $s\simeq 1/3$ have been
calibrated using measurements of velocity profiles over varied rough walls
(Schultz & Flack, 2009; Flack & Schultz, 2010) (figure 2b). In the smooth
regime ($\mathcal{R}_{d}\ll 1$), these numbers give a hydrodynamic roughness
controlled by the viscous length: $z_{0}\simeq\nu/(7u_{*})$. In the rough
regime ($\mathcal{R}_{d}\gg 1$), it rather gives $z_{0}\simeq d/30$.

The model needs to account for the hydrodynamic response of the flow to the
perturbation generated by the surface elevation profile, in the van Driest
exponential term. Let us call $\lambda$ the typical horizontal length scale of
this profile – $\lambda$ will later be the wavelength of a particular Fourier
mode. When $\lambda$ is small compared to the viscous sub-layer thickness, the
flow perturbation induced by this topography stays confined in the viscous
layer. In contrast, when $\lambda$ is much larger than $\nu/u_{*}$, the
perturbation penetrates the outer turbulent zone so that the flow response is
turbulent. In the transitional regime between laminar and turbulent responses,
the surface bumps generate turbulent perturbations that develop downstream,
due to an adverse pressure gradient (diverging mean streamlines). This
effectively modulates the thickness of the viscous sub-layer. In the mixing
length expression (5), this translates into a non-constant value of
$\mathcal{R}_{t}$. To describe these variations, Hanratty (1981) proposed an
empirical relaxation equation in which the viscous sub-layer thickness lags
behind the pressure gradient by a distance scaling with the viscous length.
Here we write it in terms of $\mathcal{R}_{t}$ as

$\frac{u_{*}}{\nu}(\mathcal{R}_{t}-\mathcal{R}_{t}^{0})+a\partial_{x}\mathcal{R}_{t}=b\frac{\mathcal{R}_{t}^{0}}{\rho
u_{*}^{2}}\partial_{x}(\tau_{xx}-p)\,.$ (6)

Both constants $a\simeq 2000$ and $b\simeq 35$ have been calibrated in Charru
et al. (2013) on measurements of the basal shear stress over a modulated
surface (Zilker et al., 1977; Zilker & Hanratty, 1979; Frederick & Hanratty,
1988). The predicted amplitude of the flow response to topography is very
sensitive to these values. In the smooth case, this response presents strong
variations for a narrow range of $\lambda$, at the transition from a viscous
response to a turbulent response (Claudin et al., 2017). In the rough case,
this transitional regime disappears as the surface roughness destabilises the
viscous sub-layer. We will show that this hydrodynamic ‘anomaly’ associated
with the laminar-turbulent transition, has a signature in the outer
hydrodynamic roughness.

###  2.3 Weakly non-linear expansion

We consider a solid surface of vertical elevation profile $Z(x)$, which is
decomposed in Fourier modes. The computation of the outer hydrodynamic
roughness induced by $Z(x)$ from the flow response to this topographical
perturbation is a non-linear calculation, as one needs to compute the
correction to the base (uniform) flow. We perform this calculation below,
starting from the base flow and then expanding the governing equations at
second order in amplitude of the $Z$ variations, considering only their
influence on the zero mode, i.e. on the average flow. Without loss of
generality, a single Fourier mode of wavenumber $k=2\pi/\lambda$ and amplitude
$\zeta$ can be considered:

$Z(x)=\zeta e^{ikx}.$ (7)

$\zeta$ is the amplitude of the surface modulation, and the following
perturbation theory will be developed in powers of $k\zeta$ assumed small. We
use standard notations with complex numbers for all quantities involved for
the linear development in the sake of mathematical convenience, though only
real parts are understood. We introduce the dimensionless vertical coordinate
$\eta=kz$, the wavenumber-based Reynolds number $\mathcal{R}={u_{*}}/{k\nu}$,
the rescaled equivalent grain diameter
$\eta_{d}=kd=\mathcal{R}_{d}/\mathcal{R}$, and the dimensionless mixing length
$\Upsilon=k\ell$.

The base state corresponds to a homogenous substrate in the $x$-direction
($\zeta=0$). In this case, the strain rate reduces into $\partial_{z}u_{x}$,
and correspondingly Eq. 3 reduces to $\partial_{z}\tau_{xz}=0$, or
equivalently $\tau_{xz}\equiv\rho|u_{*}|u_{*}$ once integrated. The mixing
length then simplifies into:

$\Upsilon=\kappa(\eta+r\eta_{d})\left[1-\exp\left(-\frac{\mathcal{R}(\eta+s\eta_{d})}{\mathcal{R}_{t}^{0}}\right)\right].$
(8)

Summing up the turbulent and viscous contributions of the shear stress
$\tau_{xz}$, (4) yields:

$\rho\left(\ell^{2}|\partial_{z}u_{x}|\partial_{z}u_{x}+\nu\partial_{z}u_{x}\right)=\rho|u_{*}|u_{*}\,.$
(9)

We define the function $\mathcal{U}(\eta)$ giving the flow velocity profile in
the base state as $u_{x}\equiv u_{*}\mathcal{U}$. Following the above
equation, $\mathcal{U}$ obeys:

$\Upsilon^{2}|\mathcal{U}^{\prime}|\mathcal{U}^{\prime}+{\mathcal{R}}^{-1}\mathcal{U}^{\prime}=1,\qquad{\rm
or\
equivalently}\qquad\mathcal{U}^{\prime}=\frac{-1+\sqrt{1+4\Upsilon^{2}{\mathcal{R}}^{2}}}{2\Upsilon^{2}{\mathcal{R}}}\,,$
(10)

where primed quantities denote derivatives with respect to $\eta$. This
equation is solved with the boundary condition $\mathcal{U}(0)=0$
corresponding to the no-slip condition of the fluid at the solid surface. The
viscous sublayer is described by $\mathcal{U}^{\prime}=\mathcal{R}$ in the
limit of negligible $\Upsilon^{2}{\mathcal{R}}^{2}$. Figure 2b shows a typical
velocity profile of a base state, with the characteristic outer logarithmic
dependence in $z$ from which the hydrodynamic roughness $z_{0}$ is extracted.
Data available in the literature are collected to demonstrate that the
predicted velocity profiles by the present model agree well with experimental
measurements.

Next, we seek to compute the flow in response to the surface perturbation (7)
one step further than in previous papers by Jackson & Hunt (1975), Sykes
(1980), Hunt et al. (1988), Fourrière et al. (2010) and Claudin et al. (2017).
At the linear order, Fourier modes can be treated independently of each other,
and the computation of the response to a pure sinusoidal profile (7) provides
the complete field. At the next quadratic order in $(k\zeta)^{2}$, the
computation rules when taking products for non-linear contributions show that
these quadratic terms are of two kinds, those in $e^{2ikx}$ and those
homogeneous in $x$. Here, only the latter, which provide corrections to the
base flow, matters. For our purpose, we consequently write the weakly
nonlinear expansion for $u_{x}$ as

$u_{x}=u_{*}\left[{\mathcal{U}}+(k\zeta)e^{ikx}U_{1}+(k\zeta)^{2}U_{0}\right]\,,$
(11)

where $U_{1}$ and $U_{0}$ are the modal functions respectively for the linear
and homogeneous quadratic responses. The systematic expansion of Navier-Stokes
equations and of the Reynolds averaged closure around the undisturbed profile
${\mathcal{U}}$ is summarised in Appendix A. An important point, beside
technical details, is that the validity of this expansion has a limited range
of $k\zeta$, which is not known in advance. Within the same expansion
framework, different representations of the solution can be built, which have
different ranges of validity (see Appendix B). In general, the elevation
profile $Z$ presents different Fourier modes, and their contributions at the
linear and homogeneous quadratic orders are additive. Importantly, even if all
modes contribute to build the hydrodynamic roughness, they do not interact. To
account for the small-scale surface corrugations one could have for example
considered them as part of the elevation profile, and worked with two well-
separated modes. Here, however, the surface roughness is encoded in the
turbulent mixing length (5) by means of two terms expressed with the
equivalent sand grain size $d$, and it affects the base velocity profile
${\mathcal{U}}$, the turbulent mixing at the surface and the homogeneous flow
response $U_{0}$. One could then wonder whether these two ways are equivalent,
and we shall see in the next section that mode coupling is required to obtain
self-consistent results and to recover the smooth-rough transition.

Figure 3: (a) Vertical profiles of the homogeneous second order correction to
the horizontal velocity for different values of $\mathcal{R}^{-1}$ and
$\mathcal{R}_{d}$ (legend), where $-U_{0}$ is shown as a function of
$\eta=kz$. (b) Homogeneous part of the velocity
$\mathcal{U}+(k\zeta)^{2}\tilde{U}_{0}$ as a function of the distance to the
surface $k(z-Z)$ for three values of $k\zeta$ (legend) and for
$\mathcal{R}^{-1}=0.05$ and $\mathcal{R}_{d}=10$. The extrapolation to $0$ of
the upper part of the profiles gives the values of the effective hydrodynamic
roughness $kz_{e}$ seen far from the surface.

###  2.4 Self-consistent looped expansion

A simple way of coupling the different modes of the elevation profile is to
ensure that the base flow used for the expansion incorporates the corrections
due to the corrugations of the solid surface at all scales. This is
fundamentally different from the introduction of new terms in the mixing
length. In practice, we use the following iterative procedure. For each mode
of given $k\zeta$ and $\mathcal{R}$, we compute the quadratic correction
$U_{0}$ to the flow induced by this perturbation using the base flow velocity
profile $\mathcal{U}$ as described in the previous paragraph. We deduce the
corrected homogeneous velocity profile
$\hat{\mathcal{U}}=\mathcal{U}+(k\zeta)^{2}U_{0}$ and compute the next
correction using $\hat{\mathcal{U}}$ (instead of $\mathcal{U}$) for the new
base flow in the weakly non-linear expansion described in Appendix A: all
convective turbulent mixing terms are then computed with a velocity profile
that accounts for that perturbation. We repeat this loop until convergence to
a final profile $\hat{\mathcal{U}}$ which, used in the expansion, leads to the
profile of $U_{0}$ giving that $\hat{\mathcal{U}}$. This self-consistent
expansion around the disturbed average velocity profile, is a priori more
precise than the non-looped expansion. This approach resembles, in spirit,
implicit (as opposed to explicit) schemes in numerical solvers of differential
equations. The price to pay, however, is this iterative procedure that must be
performed for _each_ value of $k\zeta$.

##  3 Results

We present in this section the results obtained with the above model to
compute the effective hydrodynamic roughness $z_{e}$ induced by the surface
elevation profile (7). We start with the standard perturbative calculation,
which will allow us to discuss different regimes and to evidence the signature
of the laminar-turbulent transition on $z_{e}$. With the looped calculation,
we will then show how to match inner and outer roughnesses, with application
to the smooth-rough transition as well as to superimposed bedforms.

###  3.1 Hydrodynamic roughness from homogeneous quadratic velocity correction

Following (11), the homogeneous part of the velocity profile shows that the
base state $\mathcal{U}$ is corrected by the quadratic term
$(k\zeta)^{2}U_{0}$. The modal function $U_{0}$ associated with this
correction is illustrated in figure 3a, for various values of $\mathcal{R}$
and $\mathcal{R}_{d}$. Its vertical profile tends towards a constant
$U_{0}\sim-E$ far from the surface, approximately for $z\gtrsim\lambda$. We
can define the outer hydrodynamic roughness from this asymptotic behaviour of
the function $U_{0}$, and $E$ is accordingly referred to as _roughness
coefficient_. Using the fact that the base profile
$\mathcal{U}\sim\frac{1}{\kappa}\ln\frac{z}{z_{0}}$ in that limit of large
$z$, and identifying the corresponding behaviour of the corrected velocity
profile $\mathcal{U}+(k\zeta)^{2}U_{0}$ with a similar log-profile
$\frac{1}{\kappa}\ln\frac{z}{z_{e}}$ of effective roughness $z_{e}$, we simply
obtain

$\ln\frac{z_{e}}{z_{0}}=\kappa(k\zeta)^{2}E.$ (12)

Figure 3b shows the velocity profiles in the shifted representation (Appendix
B). One observes that $z_{e}$ indeed increases with $k\zeta$. Interestingly,
the ratio $z_{e}/z_{0}$ does not only depend on the amplitude $\zeta$ of the
surface elevation, but also on $\mathcal{R}$ and $\mathcal{R}_{d}$. We
describe and analyse below the variations of the coefficient $E$, which
encodes the dependence with respect to these two quantities.

###  3.2 Roughness coefficient $E$ in the different regimes

The variations of the coefficient $E$ on the parameters $\mathcal{R}$ and
$\mathcal{R}_{d}$ are illustrated in figure 4. For a given $\mathcal{R}_{d}$,
$E$ essentially decreases with $k\nu/u_{*}=\mathcal{R}^{-1}$ (panel a), with
$E$ systematically smaller for larger $\mathcal{R}_{d}$. Below
$k\nu/u_{*}\simeq 10^{-2}$, the variations are rather weak (but note the
vertical log-scale). Above $k\nu/u_{*}\simeq 10^{-1}$, the curves become
straight and steeper, corresponding to $E\propto\mathcal{R}$, with a collapse
at small $\mathcal{R}_{d}$. When plotted as a function of $kz_{0}$, the
coefficient $E$ is also generally decreasing, except within an intermediate
range $kz_{0}\simeq 10^{-4}$–$10^{-1}$, where non-monotonic variations
develop, all the larger for smaller $\mathcal{R}_{d}$ (panel b). Below
$kz_{0}\simeq 10^{-4}$, a collapse along a parabolic behaviour is observed.
This non-monotonic behaviour is in fact also noticeable in panel (a), but less
visible due to the vertical log-scale. Below, we comment in more detail these
different regimes, and provide, when possible, some asymptotic scaling
behaviours.

Figure 4:  (a) Roughness coefficient $E$ as a function of
$k\nu/u_{*}=\mathcal{R}^{-1}$ for different $\mathcal{R}_{d}$ (legend). A line
collapse is observed in the limit of large wavenumber when $\mathcal{R}_{d}$
is small, which is well fitted with $E\sim\mathcal{R}$ (dash-dotted orange
line). Dashed lines correspond to the moment at which $\lambda$ becomes
smaller than $d$. Inset: $E/\mathcal{R}$ as a function of $\mathcal{R}_{d}$ in
the viscous limit ($\mathcal{R}\to 0$). It is well fitted with
$E/\mathcal{R}\sim\mathcal{R}_{d}^{-1}$ when $\mathcal{R}_{d}\gtrsim 100$
(dash-dotted orange line). (b) The same coefficient $E$ but now as a function
of $kz_{0}$. A line collapse is observed in the limit of small wavenumber. The
dash-dotted orange curve is the fit of the rough curves by Eq. 15, with
$a\simeq 2.7$, $b\simeq-4.1$ and $c\simeq 0$.

In the viscous regime, the flow perturbation is imbedded in the viscous sub-
layer. In this limit $\mathcal{R}$ is small and one can see in figure 4a that
this regime starts at $k\nu/u_{*}\gtrsim 10^{0}$, which roughly corresponds to
$kz_{0}\gtrsim 10^{-1}$. Since turbulent fluctuations are not relevant in this
limit, the mixing length vanishes and the shear stress is dominated by the
viscous term $\tau_{xz}\sim\rho\nu\partial_{z}u_{x}$. It follows that all
stress modal functions are simply related to the derivatives of those of the
velocity as: $S_{t1}\sim\mathcal{R}^{-1}U_{1}^{\prime}$ and
$S_{t0}\sim\mathcal{R}^{-1}U_{0}^{\prime}$. The non-linear correction to the
stress $S_{t0}$ can only come from the inertial term in the momentum balance
of Navier-Stokes equations, which becomes negligible in this limit. As a
consequence, $U_{0}^{\prime}\to 0$, so that $U_{0}$ tends to a constant, whose
value can be deduced from the lower boundary condition. In order to get no
slip at the surface, i.e. $u_{x}(Z)=0$ at second order in $k\zeta$, $U_{0}$
must compensate the variation of $U_{1}$: $-U_{0}(0)\sim U_{1}^{\prime}(0)$.
Relating this derivative to the stress, we then obtain
$-U_{0}(0)\sim\mathcal{R}S_{t1}(0)$. A proper derivation of the equations in
this asymptotic regime (see Appendix C) leads to

$E=-U_{0}\sim\frac{1}{2}\mathcal{R}\mathcal{A}\,,$ (13)

where the basal shear coefficients $S_{t1}(0)=\mathcal{A}+i\mathcal{B}$ have
been introduced (Fourrière et al., 2010; Claudin et al., 2017). This
asymptotic scaling with respect to $\mathcal{R}$ is nicely verified by the
full model (orange line in figure 4a). Furthermore, all curves collapse in
this viscous limit for the smaller values of $\mathcal{R}_{d}$. This is
consistent with the fact that the in-phase stress coefficient $\mathcal{A}\to
2$ in the smooth regime (Benjamin, 1959; Charru & Hinch, 2000; Charru et al.,
2013), so that one expects that $E\approx\mathcal{R}$, which is indeed what we
obtain up to $\mathcal{R}_{d}\simeq 100$.

For larger values of $\mathcal{R}_{d}$, the curves in figure 4a at large
wavenumber decollapse but remains parallel in this log-log representation.
This means that they share the same dependence $\propto\mathcal{R}$, but with
a prefactor that decreases with $\mathcal{R}_{d}$. Having both small
$\mathcal{R}$ and large $\mathcal{R}_{d}$ may appear to be an unphysical
limit, as it would correspond to a surface perturbation at a scale $\lambda$
smaller than the equivalent roughness size $d$. However, one could for example
imagine non-geometrical sources of effective surface roughness, such as random
jet injection (Park & Choi, 1999; Kametani et al., 2015), that could obey this
scale hierarchy. In this case, a simple argument is to keep the scaling
$E\sim\mathcal{R}$, but replacing the viscosity involved in the definition of
$\mathcal{R}$ by an effective turbulent viscosity $\nu_{t}\sim u_{*}d$, such
that

$E\sim\frac{u_{*}}{k\nu_{t}}\sim\mathcal{R}/{\mathcal{R}_{d}}\,.$ (14)

This decreasing behaviour in $1/\mathcal{R}_{d}$ at large $\mathcal{R}_{d}$ is
nicely verified by the dashed part of the curves in figure 4a (inset).

Let us turn now to the fully turbulent regime, typically for
$k\nu/u_{*}\lesssim 10^{-3}$, corresponding to $kz_{0}\lesssim 10^{-4}$. In
this limit, one can infer an expression for $E(kz_{0})$ estimating $U_{0}$ at
a typical (dimensionless) height $\eta\sim 1$, above which its profile becomes
constant (figure3a). For that purpose, we shall combine estimations of the
shear stress from the momentum balance (3) and from the turbulent closure (4).

The momentum balance essentially writes $\rho
u\partial_{x}u\sim\partial_{z}\tau$. Because the quadratic response we are
interested in is homogeneous in $x$, the only way to contribute to the
homogeneous stress correction is $U_{1}^{2}\sim S_{t0}$, where we have also
used the fact that horizontal and vertical derivatives both scale as
$\partial_{x}\sim\partial_{z}\sim 1/\eta\sim 1$. Similarly, the relation
between stress and strain rate (4) can be expressed in a scaling way as
$\tau\sim\rho\ell^{2}\left(\nabla u\right)^{2}$. For the homogeneous
correction in $(k\zeta)^{2}$, the LHS of this relation is again $S_{t0}$, but
we expect for its RHS several contributions from the product of the different
factors, combining their expressions in the base state, at linear order $1$
and quadratic homogeneous order $0$. Gathering these contributions, and as
more detailed in Appendix C, one essentially obtains $S_{t0}\sim
U_{0}+U_{1}^{2}+U_{1}+\mbox{Cst}$. Finally, using the relation for the
asymptotic expression for the linear correction of the velocity (Fourrière et
al., 2010)
$U_{1}\sim(\mathcal{A}+i\mathcal{B})/(2\kappa)\ln(\eta/kz_{0})\sim\ln(kz_{0})$,
we finally obtain for $E\sim-U_{0}$ at $\eta\sim 1$:

$\displaystyle E\simeq a\ln^{2}kz_{0}+b\ln kz_{0}+c.$ (15)

The coefficients $a$, $b$ and $c$ are expected to be approximately constant,
since they come from $\mathcal{A}$ and $\mathcal{B}$, which are only weakly
dependent on $kz_{0}$ in the turbulent regime (Fourrière et al., 2010; Charru
et al., 2013). Their adjustment to the result of the full integrated model
gives a good fit to the data for hydrodynamically rough conditions (figure
4b), in practice valid on the entire range of $kz_{0}$. Although not directly
expressed in terms of such a coefficient $E$, a similar quadratic behaviour
was derived by Jacobs (1989) by means of asymptotic analysis of RANS equations
closed by Launder-Spalding turbulence model. Taylor et al. (1989) also report
linear variations of the equivalent of $E$ with $\ln kz_{0}$ in the range
$10^{-8}$–$10^{-3}$, with different turbulent closures.

Finally, at intermediate wavenumbers, one observes a non-monotonic variation
of $E$ with respect to $kz_{0}$. It significantly increases the effective
roughness around $kz_{0}\simeq 10^{-3}$. It is most pronounced when the flow
is hydrodynamically smooth and it becomes less and less significant as
$\mathcal{R}_{d}$ increases (figure 4b). This striking behaviour is similar to
what is also observed in the flow linear response for the behaviour of the
stress coefficients $\mathcal{A}$ and $\mathcal{B}$ (Charru et al., 2013;
Claudin et al., 2017), and the pressure coefficients $\mathcal{C}$ and
$\mathcal{D}$ (Claudin et al., 2021), when $10^{-4}\lesssim kz_{0}\lesssim
10^{-2}$. This ‘anomaly’ takes place at the transition between a laminar and a
turbulent response of the flow to the elevation profile, discussed above. It
is included in the description by the relaxation equation of the transitional
Reynolds number $\mathcal{R}_{t}$ (Eq. 6), associated with this laminar-
turbulent transition. A qualitative picture of this transition is that
turbulent bursts are produced at the crest and develop on the downstream side
(Charru et al., 2013), when the disturbance reaches the size of the viscous
boundary layer. This effect shares some similarities with the drag crisis on a
sphere or a circular cylinder (Choi et al., 2008). A deeper understanding of
it is however missing and new experiments as well as numerical simulations are
definitively needed to go further on this fundamental question.

###  3.3 Matching between inner and outer scales

So far, we have modelled the inner roughness by terms in the mixing length
proportional to the equivalent grain size $d$. Consider again the case of
grains at the surface of sand ripples superimposed to a large dune. At the
largest scale, the sand ripples are the source of dune-surface roughness. They
then could be also treated by a modified mixing length with a corresponding
equivalent grain size. However, one can zoom on these sand ripples and
consider them as geometrical disturbances to the dune surface, rather than a
source of turbulent mixing. Similarly, zooming on the grains, one expects the
roughness they induce to be an emergent property of their geometrical effect.
The model should therefore be self-consistent when considering the induced
roughness from the two perspectives: a source of turbulent mixing and a
geometrical corrugation.

As a first test, we consider a flat rough surface whose corrugations are
modelled by a smooth surface of wavelength $\lambda$, which will be set by the
equivalent grain size $d$ (figure 5b). However, to avoid confusion, we will
keep the notation $\lambda$. Following the perturbative prediction (12) for a
smooth surface, the effective hydrodynamic roughness induced by such an
elevation of amplitude $\zeta$ reads

$z_{e}\simeq\frac{\nu}{7u_{*}}\exp\left[\kappa(k\zeta)^{2}E\right],$ (16)

where $E$ is given by the red curve in figure 4a (smooth limit), computed for
$k=2\pi/d$. When plotted as a function of $du_{*}/\nu$ to compare with the
data (dashed line in figure 2a, for $k\zeta=0.88$, see below), we see that the
agreement is good in the smooth regime, but starts to diverge from the
measurements when the roughness becomes on the order of the viscous sublayer:
$z_{e}u_{*}/\nu\simeq 1$. This expansion is thus unable to reproduce the
transition towards the rough regime, and predicts a hydrodynamic roughness
significantly larger than that observed: the dependence is exponential rather
than linear. The problem lies in the fact that, when the surface corrugation
increases, it induces inertial mixing of momentum that destroys the viscous
sublayer. The smooth base profile with its viscous sublayer is thus not the
relevant uniform state around which the expansion should be performed.

The idea is then to use instead the looped expansion introduced in section
2.4. We start with the the smooth base velocity profile $\mathcal{U}^{0}$,
i.e. computed from (Eq. 10) with $\eta_{d}=0$ in (8), choose a value of
$k\zeta$, and iterate the described procedure until convergence. To avoid
numerical instabilities, we use the numerical trick to approximate at each
loop $\hat{\mathcal{U}}$ by the base profile $\mathcal{U}^{d_{e}}$ computed
for a non-vanishing $d_{e}$, so that their hydrodynamic roughnesses are the
same (figure 5a).

Figure 5: (a) Velocity profiles $\mathcal{U}^{0}$, $\hat{\mathcal{U}}$ and
$\mathcal{U}^{d_{e}}$ introduced in the self-consistent calculation. These two
last functions overlap in the upper part of their profiles, i.e. possess the
same hydrodynamic roughness $z_{e}$ (dashed orange lines). Here computation
done for $k\zeta=0.88$ (see figure 6a) and $k\nu/u_{*}=0.05$, which converges
to $d_{e}u_{*}/\nu\simeq 130$. (b) Sinusoidal approximation of the surface of
a periodic array of touching grains, setting $\lambda=d$. (c) Isocontours of
the dimensionless stream function $k\Psi/u_{*}$ (see Appendix B ) above a
modulated surface (in white), when $k\zeta=0.88$ and $k\nu/u_{*}=0.05$. The
flow is from left to right. Color code: $k\Psi/u_{*}$ is increasing from blue
to red in the background; and the dashed orange curve corresponds to
$k\Psi/u_{*}=0$. Figure 6:  (a) Dimensionless effective hydrodynamic roughness
$z_{e}u_{*}/\nu$ as a function of $\lambda u_{*}/\nu$, computed in a self-
consistent way for different values of $k\zeta$ (see legend). (b) Resulting
equivalent grain size $d_{e}u_{*}/\nu$ as a function of $k\nu/u_{*}$ for
different $k\zeta$ (same colour code as in panel a).

We show in figure 6a the effective roughness computed in this self-consistent
manner for various values of $k\zeta$, as a function of $\lambda u_{*}/\nu$.
It can be seen that the transition from smooth to rough regimes is indeed
recovered: the roughness changes from $z_{e}u_{*}/\nu\simeq 1/7$ to a linear
increase with $\lambda$. This asymptotic increase is steeper for larger values
of $k\zeta$. When compared to the data on rough flat surfaces by setting
$\lambda=d$, this prediction fits well the measurements when $k\zeta\simeq
0.88$ (see the red solid curve in figure 2a). Interestingly, this value
corresponds to the sinusoidal approximation of the surface of a periodic array
of touching grains $Z\simeq 0.14\,d\cos\left(2\pi x/d\right)$, which gives
$k\zeta\simeq 0.14\times 2\pi\simeq 0.88$ (red line in figure 5b).

Although the direct and self-consistent derivations of the effective
hydrodynamic roughness (dash-dotted and solid red curves in figure 2a,
respectively) are both expansions at the second order in aspect ratio
$k\zeta$, their range of validity is rather different. The better performance
of the self-consistent approach is due to the fact that it takes into account
the topography-induced velocity reduction at the scale of $\lambda$.
Consistently, we can interpret the observation that the self-consistent
derivation predicts a lower $z_{e}$ than the direct calculation, as accounting
for flow recirculations at the scale of the topography reduces the effective
corrugation of the bottom (figure 5c).

Figure 7: Self-consistent calculation of the effective roughness _vs_ bedform
wavelength (a) and the equivalent grain size _vs_ bedform wavelength (b).
These are computed for a perturbation of various surface roughness
$\mathcal{R}_{d}$ encoded by different colours (see legend for the colour code
in both panels) and with an aspect ratio $k\zeta\simeq 0.25$, typically
representative of sedimentary dunes. Dashed lines correspond to the moment at
which $\lambda$ becomes smaller than $d$.

###  3.4 Superimposed sedimentary bedforms

The other situation where inner and outer roughnesses must match is that of
superimposed beforms: the outer roughness of a given level in the scale
hierarchy is the inner roughness of the next one, as we illustrated with
grains, ripples and dunes (figure 1a). For this purpose, the looped
calculation described in section 2.4 must be initiated at some finite value of
$\mathcal{R}_{d}$. The resulting curves are displayed in figure 7 for
$k\zeta=0.25$, a typical value for dunes at equilibrium (aspect ratio on the
order of $1/12$). It shows consistent increasing trends for larger and rougher
bed. The curves start to diverge from the smooth case when $\mathcal{R}_{d}$
is larger than a few units. In addition, the resulting effective roughness can
be further transformed into the Darcy friction factor (or skin friction) and
the hydraulic resistance (or form drag) of bedforms (Engelund, 1977; Zanke et
al., 2022). An empirical formula for the hydrodynamic roughness of sedimentary
bedforms has been proposed by Van Rijn (1984), which reflects the growth of
the hydrodynamic roughness from a value selected by the grain diameter $d$ to
a value selected by the bedform wavelength, as the aspect ratio $k\zeta$
increases. At small aspect ratio, for $k\zeta\lesssim 0.25$, which includes
the typical aspect ratio of sedimentary ripples and dunes, the looped
calculation performed here turns out to be recovering values close to this
empirical prediction.

This situation of multi-scaled pattern is particularly relevant in the sub-
aquatic case, as well as for low-pressure – e.g. planetary – gaseous
environments, which provide viscous or transitional conditions (Colombini &
Stocchino, 2011; Jia et al., 2017; Durán Vinent et al., 2019; Andreotti et
al., 2021; Gunn & Jerolmack, 2022). We shall illustrate here how to apply the
present results on the emblematic example of Martian dunes and ripples (figure
1b,c). The dunes, whose typical length is on the order of a few hundreds of
meters, exhibit large meter-scale ripples on their surface; these ripples
themselves are superimposed with decimeter-scale small ripples (Lapôtre et
al., 2016, 2018). A hydrodynamic-based explanation for the existence of these
intermediate bedforms has been proposed by Durán Vinent et al. (2019), and
verified by Rubanenko et al. (2022): the anomaly associated with the laminar-
turbulent transition creates a forbidden gap in wavelength that stops the
coarsening of the most unstable mode. As discussed above, this anomaly is
present in the smooth regime only. Estimating Martian conditions with a grain
size $d\simeq 100~{}\mu$m, a shear velocity $u_{*}\simeq 0.5$ m/s and a
viscosity $\nu\simeq 10^{-3}$ m2/s, one obtains $\mathcal{R}_{d}\simeq 0.05$,
i.e. a hydrodynamically smooth granular bed. The question is then to evaluate,
for the large meter-scale ripples, the hydrodynamic roughness induced by their
superimposed small ripples. Taking, for these small ripples, a wavelength
$\lambda\simeq 0.1$ m and an aspect ratio on the order of $1/20$, i.e.
$k\zeta\simeq 0.15$, one obtains, with $\lambda u_{*}/\nu\simeq 50$, an
effective roughness $z_{e}u_{*}/\nu\simeq 0.2$ corresponding to an equivalent
sand grain Reynolds number $\mathcal{R}_{d_{e}}\simeq 5$. This indicates that
the flow over the large ripples is indeed in the smooth regime, consistently
with the hypothesis of Durán Vinent et al. (2019). Now taking for these large
ripples $\lambda\simeq 3$ m, i.e. $\lambda u_{*}/\nu\simeq 1.5\times 10^{3}$,
and $k\zeta\simeq 0.25$, we see in figure 7 ($\mathcal{R}_{d}=5$ is between
the blue and black lines) an effective roughness $z_{e}u_{*}/\nu\simeq 1.5$,
corresponding to an equivalent sand grain Reynolds number
$\mathcal{R}_{d_{e}}\simeq 55$. This means that the flow over the large dunes
should be considered intermediate between smooth and rough regimes.
Importantly, aeolian saltation – and sediment transport in general (Van Rijn,
1984) – is known to increase the sediment bed roughness, with $z_{0}$ on the
order of a few $d$ in the terrestrial case, associated with Bagnold’s focal
point (Durán et al., 2011; Valance et al., 2015). In low-pressure (Martian)
conditions, no such measurement is available so far, but we have checked that
the above conclusions are unchanged when doing this hierarchical looped
calculation with an initial similar saltation-induced surface roughness.

###  3.5 Non-linear size selection of dissolution bedforms

Finally, the self-consistent prediction of the effective roughness induced by
a wavy surface can be used to revisit the argument given by Claudin et al.
(2017) for the selection of the aspect ratio of dissolution patterns. The
instability at the origin of these bedforms is associated with the
hydrodynamic anomaly on smooth modulated beds for a narrow range of
wavenumbers around $k\nu/u_{*}\simeq 10^{-3}$. This anomaly, and thus the
instability mechanism, disappears when the bed is rough enough. Because the
development of the bedforms generates an effective roughness, the idea is that
they would stop growing when this threshold in $\mathcal{R}_{d}$ is reached.
Here we can make this argument more quantitative with the proposed non-linear
calculations. The looped computation of $z_{e}$ equivalently provides a
corresponding equivalent grain size $d_{e}$, which depends on $k\zeta$ (figure
6b). Experimental studies on the development of ice ripples and scallops
report equilibrium shapes with the inverse Reynolds numbers in the range
$k\nu/u_{*}\simeq 8\times 10^{-4}$ – $2\times 10^{-3}$, and aspect ratios
$6$–$8$ %, i.e. $k\zeta\simeq 0.19$–$0.25$ (Ashton & Kennedy, 1972; Bushuk et
al., 2019). The corresponding induced roughness reads from figure 6b around
$d_{e}u_{*}/\nu\simeq 50$. This value is nicely consistent with the
instability threshold computed for parameters relevant to the case of ice
melting (Claudin et al., 2017).

##  4 Concluding remarks

Pursuing our goal to describe the hydrodynamic response to a bed perturbation,
we have here gone one step beyond the linear order, and computed the
homogeneous correction to the base flow induced by the topography. This is a
quadratic effect in the amplitude of the bed elevation, thus requiring weakly
non-linear calculations. Three main results can be emphasised: we quantify how
the bed corrugation increases the effective hydrodynamic roughness; we show
that this effect is sensitive to the laminar-turbulent transition; this model
is able to reproduce the smooth-rough transition in a self-consistent manner
and can be applied to multiscale elevation profiles where inner and outer
roughnesses are hierarchically nested.

The expression (12) shows that the effective roughness $z_{e}$ cannot be
simply related to a single geometrical length (Schlichting & Gersten, 2000;
Van Rijn, 1982; Raupach et al., 1991a; Wiberg & Nelson, 1992), here the
amplitude $\zeta$ of the surface elevation profile. It still depends on the
inner roughness of the flat surface $z_{0}$. Moreover, it involves a
coefficient $E$ that also depends on the wavenumber of the surface profile,
typically decreasing with larger $k$. Interestingly, it shows a non-monotonic
behaviour for a range of wavenumber. This anomaly is larger for smoother
surfaces, and disappears in the rough limit. Consistently with a similar
anomaly for the stress and pressure response (Claudin et al., 2017, 2021), it
can be associated, in the model, to the lag between the thickness of the
viscous sub-layer (parametrised by $\mathcal{R}_{t}$) and the pressure
gradient (Eq. 6). This anomaly plays a key role in the formation of
sedimentary, dissolution or sublimation bedforms, especially in low-pressure
planetary conditions (Durán Vinent et al., 2019; Bordiec et al., 2020).
However, its direct experimental evidence is limited to a series of
measurements by Hanratty and co-workers (Zilker et al., 1977; Frederick &
Hanratty, 1988), that involve a dedicated apparatus with electro-kinetic
probes to obtain the basal shear stress response. The present work suggests
that the measurement of the hydrodynamic roughness induced by a sinusoidal
smooth bottom, achievable by more standard and non-intrusive velocimetric
techniques, could offer an alternative way to provide experimental evidence
for this anomaly and allow for a better calibration of the relaxation equation
(6). Importantly, as a larger surface perturbations induce the smooth to rough
transition, the self-consistent calculation of $z_{e}$ shows that the anomaly
disappears when $k\zeta$ is too large, and the non-monotonic behaviour of the
roughness is effectively observed here typically for $k\zeta\lesssim 0.05$.
Finally, these results also shed light on open hydrodynamic topics associated
with such transitional shear flows over a solid wall (Tuckerman et al., 2020;
Gomé et al., 2022).

Numerical simulations provide another way to investigate these questions. In
particular, DNS over wavy surfaces, in the spirit of those of Maaß & Schumann
(1996) and De Angelis et al. (1997), would allow one to gain a deeper
understanding of Hanratty’s anomaly. In such simulations, the good control of
the imposed flow as well as the surface properties (wavelength, amplitude,
roughness), together with a scale separation which requires a resolved viscous
sub-layer much smaller than the surface wavelength, itself sufficiently
smaller than the system size, is definitively challenging at relevant values
of the Reynolds number, but probably reachable with current numerical
techniques and computer power (Lee & Moser, 2015). The understanding and the
description of the interplay between a wavy surface and the associated
modulation of the viscous sublayer then remains an interesting open problem
both from numerical and experimental points of view, with significant
potential applications for geomorphology and geophysical flows.

##  Appendix A Weakly non-linear expansion

In this appendix, we provide detailed information on how to perform the weakly
nonlinear expansion. The aim is to derive two sets of closed equations for the
linear and homogeneous quadratic responses, which are then solved with
boundary conditions at the bottom and top of the computing domain.

The weakly nonlinear expansion for $u_{x}$ has been given in Eq. 11 as follows

$u_{x}=u_{*}\left[{\mathcal{U}}+(k\zeta)e^{ikx}U_{1}+(k\zeta)^{2}U_{0}\right]\,,$
(17)

with $U_{1}$ and $U_{0}$ being the modal functions respectively for the linear
and homogeneous quadratic responses. Herein, standard notations with complex
numbers are used for all quantities involved for the linear development in the
sake of mathematical convenience, though only real parts are understood.
However, when computing the non-linear terms, one must go back to real
notations through the transform $f\to(f+f^{*})/2$, with $f^{*}$ being the
complex conjugate of $f$. In the same way, expansions are respectively
performed for the vertical velocity $u_{z}$, the shear stress $\tau_{xz}$, the
normal stress difference $\tau_{zz}-\tau_{xx}$, the vertical component of the
stress $p-\tau_{zz}$ and the dimensionless mixing length $k\ell$ as

$\displaystyle u_{z}$ $\displaystyle=$ $\displaystyle
u_{*}\left[(k\zeta)e^{ikx}W_{1}\right]\,,$ (18) $\displaystyle\tau_{xz}$
$\displaystyle=$ $\displaystyle\rho
u_{*}^{2}\left[1+(k\zeta)e^{ikx}S_{t1}+(k\zeta)^{2}S_{t0}\right]\,,$ (19)
$\displaystyle\tau_{zz}-\tau_{xx}$ $\displaystyle=$ $\displaystyle\rho
u_{*}^{2}\left[(k\zeta)e^{ikx}S_{d1}+(k\zeta)^{2}S_{d0}\right]\,,$ (20)
$\displaystyle p-\tau_{zz}$ $\displaystyle=$ $\displaystyle\rho
u_{*}^{2}\left[\frac{1}{3}\chi^{2}+(k\zeta)e^{ikx}S_{n1}+(k\zeta)^{2}S_{n0}\right]\,,$
(21) $\displaystyle k\ell$ $\displaystyle=$ $\displaystyle\Upsilon+k\zeta
e^{ikx}L_{1}+(k\zeta)^{2}L_{0}\,.$ (22)

Note that as there cannot be fluid flow through the solid boundary, the term
$W_{0}(\eta)=0$ is not included into the expansion.

With these notations, the strain rate modulus can be rewritten as

$\displaystyle|\dot{\gamma}|$ $\displaystyle=$
$\displaystyle\mathcal{U}^{\prime}+(k\zeta)e^{ikx}(U_{1}^{\prime}+iW_{1})+(k\zeta)^{2}\left[U_{0}^{\prime}+\frac{1}{\mathcal{U}^{\prime}}U_{1}U_{1}^{*}\right]\,,$
(23)

and following Eq. 4, the stress modal functions can be expressed as

$\displaystyle S_{t1}$ $\displaystyle=$
$\displaystyle\left(2\Upsilon^{2}\mathcal{U}^{\prime}+\mathcal{R}^{-1}\right)(U_{1}^{\prime}+iW_{1})+2L_{1}\Upsilon\mathcal{U}^{\prime
2}\,,$ (24) $\displaystyle S_{t0}$ $\displaystyle=$
$\displaystyle\left(2\Upsilon^{2}\mathcal{U}^{\prime}+\mathcal{R}^{-1}\right)U_{0}^{\prime}+\frac{1}{2}L_{1}\mathcal{U}^{\prime
2}L_{1}^{*}+\Upsilon^{2}\left(U_{1}U_{1}^{*}+\frac{1}{2}(U_{1}^{\prime}+iW_{1})\left(U_{1}^{\prime*}-iW_{1}^{*}\right)\right)\,,$
(25) $\displaystyle+$
$\displaystyle\Upsilon\left(\mathcal{U}^{\prime}L_{1}^{*}(U_{1}^{\prime}+iW_{1})+L_{1}\mathcal{U}^{\prime}\left(U_{1}^{\prime*}-iW_{1}^{*}\right)+2L_{0}\mathcal{U}^{\prime
2}\right)\,,$ $\displaystyle S_{d1}$ $\displaystyle=$
$\displaystyle-\frac{4i}{\mathcal{U}^{\prime}}U_{1}\,,$ (26) $\displaystyle
S_{d0}$ $\displaystyle=$
$\displaystyle-2i\Upsilon\left(\mathcal{U}^{\prime}U_{1}L_{1}^{*}-L_{1}\mathcal{U}^{\prime}U_{1}^{*}\right)-i\Upsilon^{2}\left(U_{1}\left(U_{1}^{\prime*}-iW_{1}^{*}\right)-U_{1}^{*}(U_{1}^{\prime}+iW_{1})\right)\,,$
(27)

where
$(\mathcal{R}^{-1}+\Upsilon^{2}\mathcal{U}^{\prime})=1/\mathcal{U}^{\prime}$
as deduced from Eq. 10. Further, the modal functions for the mixing length are
obtained as

$\displaystyle\frac{L_{1}}{\kappa}$ $\displaystyle=$
$\displaystyle\frac{1}{2}e^{-(\eta+s\eta_{d})\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}}\left\\{(\eta+r\eta_{d})\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}\Big{[}2\left(\eta+s\eta_{d}\right)R_{t1}+\left(\eta+s\eta_{d}\right)S_{t1}-2\Big{]}+2\right\\}-1\,,$
(28) $\displaystyle\frac{L_{0}}{\kappa}$ $\displaystyle=$
$\displaystyle-\frac{1}{16}\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}\,e^{-(\eta+s\eta_{d})\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}}\left\\{4\Big{[}(\eta+r\eta_{d})\left(\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}-2S_{t0}(\eta+s\eta_{d})\right)-2\right.$
(29) $\displaystyle+$
$\displaystyle\left(R_{t1}^{*}+R_{t1}+(S_{t1}+S_{t1}^{*})/2\right)\left((\eta+r\eta_{d})+(\eta+s\eta_{d})-(\eta+r\eta_{d})(\eta+s\eta_{d})\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}\right)\Big{]}$
$\displaystyle+$ $\displaystyle
4(\eta+r\eta_{d})(\eta+s\eta_{d})R_{t1}^{*}R_{t1}S_{t1}\left((\eta+s\eta_{d})\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}-2\right)\left((\eta+s\eta_{d})\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}-1\right)$
$\displaystyle+$
$\displaystyle\left.(\eta+r\eta_{d})(\eta+s\eta_{d})S_{t1}^{*}\left[2R_{t1}\left((\eta+s\eta_{d})\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}-1\right)+(\eta+s\eta_{d})\frac{\mathcal{R}}{\mathcal{R}_{t}^{0}}S_{t1}+S_{t1}\right]\right\\}\,,$

where

$R_{t1}=\frac{ib}{\mathcal{R}+ia}(S_{d1}+S_{n1})$ (30)

is defined from the transitional Reynolds number
$\mathcal{R}_{t}/\mathcal{R}_{t}^{0}=1-(k\zeta)e^{ikx}R_{t1}$ and deduced from
Hanratty’s relaxation relation (Eq. 6).

Introducing the above expansions into the Navier-Stokes equations (Eqs. 2 \-
3), we derive the differential equations at linear and homogeneous-quadratic
orders on the modal functions. The linear components yield:

$\displaystyle U_{1}^{\prime}$ $\displaystyle=$ $\displaystyle-
iW_{1}+\frac{S_{t1}-2\Upsilon\mathcal{U}^{\prime
2}L_{1}}{{\mathcal{R}}^{-1}+2\Upsilon^{2}\mathcal{U}^{\prime}}\,,$ (31)
$\displaystyle W_{1}^{\prime}$ $\displaystyle=$ $\displaystyle-iU_{1}\,,$ (32)
$\displaystyle S_{t1}^{\prime}$ $\displaystyle=$ $\displaystyle
i\mathcal{U}U_{1}+\mathcal{U}^{\prime}W_{1}+i(S_{n1}+S_{d1})\,,$ (33)
$\displaystyle S_{n1}^{\prime}$ $\displaystyle=$
$\displaystyle-i\mathcal{U}W_{1}+iS_{t1}\,,$ (34)

where the expression giving $U_{1}^{\prime}$ (Eq. 31) results from the shear
stress modal function $S_{t1}$ (Eq. 24). This is the result derived in Claudin
et al. (2017). For the new homogeneous-quadratic order, the differential
equations of the modal functions are:

$\displaystyle U_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{\left(2\Upsilon^{2}\mathcal{U}^{\prime}+\mathcal{R}^{-1}\right)}\Bigg{[}\Bigg{.}S_{t0}-\frac{1}{2}L_{1}\mathcal{U}^{\prime
2}L_{1}^{*}-\Upsilon^{2}\left(U_{1}U_{1}^{*}+\frac{1}{2}(U_{1}^{\prime}+iW_{1})\left(U_{1}^{\prime*}-iW_{1}^{*}\right)\right)\Bigg{.}$
(35) $\displaystyle-$
$\displaystyle\Upsilon\left(\mathcal{U}^{\prime}L_{1}^{*}(U_{1}^{\prime}+iW_{1})+L_{1}\mathcal{U}^{\prime}\left(U_{1}^{\prime*}-iW_{1}^{*}\right)+2L_{0}\mathcal{U}^{\prime
2}\right)\Bigg{]}\,,$ $\displaystyle S_{t0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{4}(W_{1}U_{1}^{\prime*}+W_{1}^{*}U_{1}^{\prime})\,,$
(36) $\displaystyle S_{n0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{1}{2}i(U_{1}W_{1}^{*}-U_{1}^{*}W_{1})\,,$ (37)

where, as in the linear order, the expression for $U_{0}^{\prime}$ (Eq. 35) is
deduced from the shear stress modal function $S_{t0}$ (Eq. 25). Note that
$S_{n0}$ does not enter the equations for $U^{\prime}_{0}$ and
$S^{\prime}_{t0}$, and it thus decouples from the present problem focused on
the computation of $U_{0}$. Furthermore, the stress functions $S_{t0}$ and
$S_{n0}$ can be obtained by integration over $\eta$, once the linear order
($U_{1}$ and $W_{1}$) is known.

The integration of the two above sets of differential equations (31)-(37)
requires boundary conditions at the substrate surface and at infinity. The
upper boundary corresponds to the limit $\eta\to\infty$, in which the vertical
flux of momentum vanishes asymptotically. This means that the corrections to
the vertical velocity and to the shear stress must tend to zero at both
orders: (i) $W_{1}(\infty)=0$ and (ii) $S_{t1}(\infty)=0$, $S_{t0}(\infty)=0$.
In practice, we introduce a finite height $H$ (or $\eta_{H}\equiv kH$), at
which we impose a null vertical velocity and a constant tangential stress
$-\rho u_{*}^{2}$, so that $W_{1}(\eta_{H})=0$ and $S_{t1}(\eta_{H})=0$,
$S_{t0}(\eta_{H})=0$. Then, we consider the limit $H\to+\infty$, i.e. when the
results become independent of $H$. Both components of the velocity should
vanish at the substrate surface, i.e. for $\eta=kZ$. The condition
$u_{x}(x,kZ)=0$ then yields:

$\displaystyle U_{1}(0)$ $\displaystyle=$
$\displaystyle-\mathcal{U}^{\prime}(0)\,,$ (38) $\displaystyle U_{0}(0)$
$\displaystyle=$
$\displaystyle-\frac{1}{4}\mathcal{U}^{\prime\prime}(0)-\frac{1}{4}U_{1}^{\prime}(0)-\frac{1}{4}U_{1}^{\prime*}(0).$
(39)

Similarly, for $u_{z}(x,kZ)=0$, we get $W_{1}(0)=0$.

The differential equations (31 \- 34) for the functions associated with the
linear terms are first solved by integrating with the corresponding boundary
conditions. The solutions are then used to solve equations (35 \- 36) for the
functions associated with the homogeneous terms.

##  Appendix B Shifted representations and stream function

A ‘shifted’ representation of the modal functions is used in some parts of the
paper, noted with an additional tilde. In this study, all fields are expanded
up to the homogeneous second order in $k\zeta$ in a non-shifted representation
as in Eq. 11 for the streamwise velocity. For some of the graphs, we sometimes
prefer to present the velocity profiles in the shifted representation as

$u_{x}={u_{*}}\left[\mathcal{U}(\xi)+(k\zeta)e^{ikx}\tilde{U}_{1}(\xi)+(k\zeta)^{2}\tilde{U}_{0}(\xi)\right]\,,$
(40)

where $\xi=\eta-k\zeta e^{ikx}$. Expanding the functions
$\tilde{f}(\eta-k\zeta e^{ikx})$ with respect to $k\zeta$, one obtains

$\tilde{U}_{0}=U_{0}+\frac{1}{4}\mathcal{U}^{\prime\prime}+\frac{1}{4}\left(U_{1}^{\prime}+U_{1}^{\prime*}\right)\,,$
(41)

where $\tilde{U}_{1}=U_{1}+\mathcal{U}^{\prime}$ has been given in Fourrière
et al. (2010). The shifted representations for the other fields work in a
similar fashion, except that there could be no base profile as $\mathcal{U}$
in the expressions.

To compute the streamlines, we introduce the stream function $\Psi(x,z)$. In
the shifted representation, one solution is $\Psi=\int_{0}^{z}u_{x}{\rm
d}\hat{z}$, and this integral is computed between the surface $\hat{z}=0$ and
some vertical distance $\hat{z}=z$. Following Eq. 40, one then has

$\Psi=\frac{u_{*}}{k}\int_{0}^{\xi}\left[\mathcal{U}(\hat{\xi})+(k\zeta)e^{ikx}\tilde{U}_{1}(\hat{\xi})+(k\zeta)^{2}\tilde{U}_{0}(\hat{\xi})\right]{\rm
d}\hat{\xi}\,.$ (42)

Introducing $\tilde{U}_{1}=U_{1}+\mathcal{U}^{\prime}$ and Eq. 41, it gives

$\Psi=\frac{u_{*}}{k}\int_{0}^{\eta}\left[\mathcal{U}+(k\zeta)^{2}\left(U_{0}+\frac{1}{4}\mathcal{U}^{\prime\prime}+\frac{1}{4}\left(U_{1}^{\prime}+U_{1}^{\prime*}\right)\right)+(k\zeta)e^{ikx}\left(iW^{\prime}_{1}+\mathcal{U}^{\prime}\right)\right]{\rm
d}\hat{\eta}\,,$ (43)

where we have used the relation at the linear order $W_{1}^{\prime}=-iU_{1}$
(Eq. 33). Consequently, one obtains a dimensionless stream function
$k\Psi/u_{*}$

$\displaystyle\frac{k\Psi}{u_{*}}$ $\displaystyle=$
$\displaystyle\int_{0}^{\eta}\left(\mathcal{U}+(k\zeta)^{2}U_{0}\right){\rm
d}\hat{\eta}+\frac{1}{4}(k\zeta)^{2}\left[\mathcal{U}^{\prime}+\left(U_{1}+U_{1}^{*}\right)+\mathcal{U}^{\prime}(0)\right]$
(44) $\displaystyle+$
$\displaystyle(k\zeta)e^{ikx}\left(iW_{1}+\mathcal{U}\right),$

where we have used the boundary conditions $\mathcal{U}(0)=0$, $W_{1}(0)=0$
and $U_{1}(0)=-\mathcal{U}^{\prime}(0)$.  For the stream function in the self
- consistent expansion analysis (section 2.4), it should be noted that in this
equation the homogeneous velocity of the base state $\mathcal{U}$ in the
integration should be estimated with $\mathcal{R}_{d}=0$ as $\mathcal{U}^{0}$,
and that the rest quantities should be estimated with
$\mathcal{R}_{d}=d_{e}u_{*}/\nu$.  Here, for the plots (figure 5c), only the
real part of (44) is considered.

##  Appendix C Asymptotic regimes of the roughness coefficient $E$

In this appendix, we provide more details in the scaling arguments for the
viscous and turbulent regimes of the roughness coefficient $E$.

###  C.1 Viscous regime

In the viscous response limit, the flow is governed by the Stokes equations.
At the linear order in $k\zeta$, the corresponding equations on the modal
functions are

$\displaystyle U_{1}^{\prime}$ $\displaystyle=$ $\displaystyle-
iW_{1}+\frac{S_{t1}}{{\mathcal{R}}^{-1}}\,,$ (45) $\displaystyle
W_{1}^{\prime}$ $\displaystyle=$ $\displaystyle-iU_{1}\,,$ (46) $\displaystyle
S_{t1}^{\prime}$ $\displaystyle=$ $\displaystyle i(S_{n1}+S_{d1})\,,$ (47)
$\displaystyle S_{n1}^{\prime}$ $\displaystyle=$ $\displaystyle iS_{t1}\,.$
(48)

At the homogeneous quadratic order, we have

$\displaystyle U_{0}^{\prime}$ $\displaystyle=$
$\displaystyle\frac{S_{t0}}{\mathcal{R}^{-1}}\,,$ (49) $\displaystyle
S_{t0}^{\prime}$ $\displaystyle=$ $\displaystyle 0\,,$ (50) $\displaystyle
S_{n0}^{\prime}$ $\displaystyle=$ $\displaystyle 0\,.$ (51)

Using the boundary condition $S_{t0}(\eta_{H})=0$, we can deduce from (50)
that $S_{t0}=0$. From (49), we therefore obtain $U^{\prime}_{0}=0$, i.e.
$U_{0}$ is a constant. Considering the no-slip boundary condition at the
surface (Eq. 39), we have

$U_{0}=U_{0}(0)=-\frac{1}{4}\mathcal{U}^{\prime\prime}(0)-\frac{1}{4}U_{1}^{\prime}(0)-\frac{1}{4}U_{1}^{\prime*}(0)\,.$
(52)

In this viscous limit $\mathcal{U}^{\prime}=\mathcal{R}$, independent of
$\eta$, and therefore $\mathcal{U}^{\prime\prime}(0)=0$ . From (45), we deduce

$U_{1}^{\prime}(0)={S_{t1}(0)}{{\mathcal{R}}}\,,$ (53)

where we have used the condition $W_{1}(0)=0$. Considering
$S_{t1}(0)=\mathcal{A}+i\mathcal{B}$, we then have

$U_{0}=-\frac{1}{4}(S_{t1}(0)+S_{t1}^{*}(0))=-\frac{1}{2}\mathcal{R}\,\mathcal{A}\,,$
(54)

and thus

$E=-U_{0}=\frac{1}{2}\mathcal{R}\,\mathcal{A}\,.$ (55)

###  C.2 Turbulent regime

In this limit, one can infer an expression for $E(kz_{0})$ estimating $U_{0}$
at a typical (dimensionless) height $\eta\sim 1$, above which its profile
becomes constant (figure3a). For that purpose, we combine estimations of the
shear stress from the momentum balance (3) and from the turbulent closure (4).

Let us start with (4), we can express in this limit as
$\tau\sim\rho\ell^{2}\left(\nabla u\right)^{2}$. For the homogeneous
correction in $(k\zeta)^{2}$, it writes

$S_{t0}=\Upsilon^{2}\mathcal{U}^{\prime}U_{0}^{\prime}+L_{0}^{2}\mathcal{U}^{\prime
2}+\Upsilon^{2}U_{1}^{2}+\Upsilon^{2}U_{1}^{\prime 2}+\Upsilon
L_{1}\mathcal{U}^{\prime}U_{1}^{\prime}+L_{1}^{2}\mathcal{U}^{\prime 2}$ (56)

where the first term on the RHS comes from the velocity gradient taken at the
quadratic order $U_{0}$, with all other factors taken in the base state
($\Upsilon$ for the dimensionless mixing length and $\mathcal{U}$ for the
dimensionless velocity). The second term is the modulated mixing length at
this order 0, and the base velocity. The other contributions come from
combinations of terms at linear order. A first possibility is to take a base
mixing length with velocity horizontal and vertical derivatives. A second
possibility is with a modulated mixing length and (vertical) velocity
derivative. A last possibility is to take a modulated mixing length with a
base velocity. One can recognise all these terms in Eq. 25, properly expressed
with complex notations and correct prefactors. We can estimate their scaling
behaviour having in mind that, in this rough turbulent limit, we have:
$\Upsilon\sim\kappa\eta$, $\mathcal{U}^{\prime}\sim 1/\Upsilon$,
$L_{1}\sim\kappa$ and $L_{0}\to 0$. Also, we estimate quantities and vertical
derivatives at the scale $\eta\sim 1$, so that
$\mathcal{U}^{\prime}\sim\mathcal{U}/\eta\sim\mathcal{U}$ (and similarly for
$U_{1}$ and $U_{0}$). It essentially leads to $S_{t0}\sim
U_{0}+U_{1}^{2}+U_{1}+\mbox{Cst}$.

[Acknowledgements]We thank F. Charru, O. Durán Vinent and M. Louge for long
term collaborations on this subject of hydrodynamic response to surface
topography. We acknowledge the contribution of A. Fourrière in setting a first
version of these non-linear hydrodynamic calculations. P.C. thanks W.
Anderson, S. Carpy, S. Courrech du Pont, L. Couston, L. Duchemin, B. Favier
and C. Narteau for discussions on Hanratty’s anomaly.

[Funding]Pan Jia was supported by Shenzhen Science and Technology Programme
(Grant No. 629 RCBS20200714114940144), and National Natural Science Foundation
of China (NSFC12102114).

[Declaration of interests]The authors report no conflict of interest.

[Data availability statement]The data that support the findings of this study
are available from the authors upon request.

[Author ORCID]Pan Jia, https://orcid.org/0000-0002-3988-4363; Bruno Andreotti,
https://orcid.org/0000-0001-8328-6232; Philippe Claudin,
https://orcid.org/0000-0001-8975-4500.

## References

  * Andreotti (2004) Andreotti, B. 2004 A two-species model of aeolian sand transport. J. Fluid Mech. 510, 47–70.
  * Andreotti et al. (2021) Andreotti, B., Claudin, P., Iversen, J.J., Merrison, J.P. & Rasmussen, K.R. 2021 A lower-than-expected saltation threshold at martian pressure and below. Proc. Natl. Acad. Sci. U.S.A. 118 (5), e2012386118.
  * Ashton & Kennedy (1972) Ashton, G.D. & Kennedy, J.F. 1972 Ripples on underside of river ice covers. J. Hydraul. Div. 98 (9), 1603–1624.
  * Aulnette et al. (2019) Aulnette, M., Rabaud, M. & Moisy, F. 2019 Wind-sustained viscous solitons. Phys. Rev. Fluids 4 (8), 084003\.
  * Ayotte et al. (1994) Ayotte, K.W., Xu, D.P. & Taylor, P.A. 1994 The impact of turbulence closure schemes on predictions of the mixed spectral finite-difference model for flow over topography. Boundary Layer Meteorol 68 (1), 1–33.
  * Bagnold (1941) Bagnold, R.A. 1941 The physics of blown sand and desert dunes. London: Methuen.
  * Benjamin (1959) Benjamin, T.B. 1959 Shearing flow over a wavy boundary. J. Fluid Mech. 6 (2), 161–205.
  * Bordiec et al. (2020) Bordiec, M., Carpy, S., Bourgeois, O., Herny, C., Massé, M., Perret, L., Claudin, P., Pochat, S. & Douté, S. 2020 Sublimation waves: Geomorphic markers of interactions between icy planetary surfaces and winds. Earth-Sci. Rev. 211, 103350.
  * Bradley (1980) Bradley, E.F. 1980 An experimental study of the profiles of wind speed, shearing stress and turbulence at the crest of a large hill. Q. J. R. Meteorol. Soc. 106 (447), 101–123.
  * Brown et al. (2008) Brown, S., Nickling, W.G. & Gillies, J.A. 2008 A wind tunnel examination of shear stress partitioning for an assortment of surface roughness distributions. J. Geophys. Res. Earth. Surf. 113 (F2).
  * Bushuk et al. (2019) Bushuk, M., Yang, X., Winton, M., Msadek, R., Harrison, M., Rosati, A. & Gudgel, R. 2019 The value of sustained ocean observations for sea ice predictions in the barents sea. J. Clim. 32 (20), 7017–7035.
  * Charru et al. (2013) Charru, F., Andreotti, B. & Claudin, P. 2013 Sand ripples and dunes. Annu. Rev. Fluid Mech. 45, 469–493.
  * Charru & Hinch (2000) Charru, F. & Hinch, E.J. 2000 ‘phase diagram’of interfacial instabilities in a two-layer couette flow and mechanism of the long-wave instability. J. Fluid Mech. 414, 195–223.
  * Cheng et al. (2007) Cheng, H., Hayden, P., Robins, A.G. & Castro, I.P. 2007 Flow over cube arrays of different packing densities. J. Wind Eng. Ind. Aerodyn. 95 (8), 715–740.
  * Choi et al. (2008) Choi, H., Jeon, W.P. & Kim, J. 2008 Control of flow over a bluff body. Annu. Rev. Fluid Mech. 40, 113–139.
  * Chung et al. (2021) Chung, D., Hutchins, N., Schultz, M.P. & Flack, K.A. 2021 Predicting the drag of rough surfaces. Annu. Rev. Fluid Mech. 53, 439–471.
  * Claudin et al. (2017) Claudin, P., Durán, O. & Andreotti, B. 2017 Dissolution instability and roughening transition. J. Fluid Mech. 832.
  * Claudin et al. (2021) Claudin, P., Louge, M. & Andreotti, B. 2021 Basal pressure variations induced by a turbulent flow over a wavy surface. Front. Phys. 9, 682564.
  * Colombini (2004) Colombini, M. 2004 Revisiting the linear theory of sand dune formation. J. Fluid Mech. 502, 1–16.
  * Colombini & Stocchino (2011) Colombini, M. & Stocchino, A. 2011 Ripple and dune formation in rivers. J. Fluid Mech. 673, 121–131.
  * Crago et al. (2012) Crago, R.D., Okello, W. & Jasinski, M.F. 2012 Equations for the drag force and aerodynamic roughness length of urban areas with random building heights. Bound. Layer Meteorol. 145 (3), 423–437.
  * De Angelis et al. (1997) De Angelis, V., Lombardi, P. & Banerjee, S. 1997 Direct numerical simulation of turbulent flow over a wavy wall. Phys. Fluids 9 (8), 2429–2442.
  * De Langre (2008) De Langre, E. 2008 Effects of wind on plants. Annu. Rev. Fluid Mech. 40, 141–168.
  * Dong et al. (2002) Dong, Z., Liu, X. & Wang, X. 2002 Aerodynamic roughness of gravel surfaces. Geomorphology 43 (1-2), 17–31.
  * Durán et al. (2011) Durán, O., Claudin, P. & Andreotti, B. 2011 On aeolian transport: Grain-scale interactions, dynamical mechanisms and scaling laws. Aeolian Res. 3 (3), 243–270.
  * Durán Vinent et al. (2019) Durán Vinent, O., Andreotti, B., Claudin, P. & Winter, C. 2019 A unified model of ripples and dunes in water and planetary environments. Nat. Geosci. 12 (5), 345–350.
  * Elbelrhiti et al. (2005) Elbelrhiti, H., Claudin, P. & Andreotti, B. 2005 Field evidence for surface-wave-induced instability of sand dunes. Nature 437 (7059), 720–723.
  * Engelund (1977) Engelund, F. 1977 Hydraulic resistance for flow over dunes. Prog. Rep 44, 19M20.
  * Fang & Sill (1992) Fang, C. & Sill, B.L. 1992 Aerodynamic roughness length: correlation with roughness elements. J. Wind Eng. Ind. Aerodyn. 41 (1-3), 449–460.
  * Finnigan (1988) Finnigan, J.J. 1988 Air flow over complex terrain. Flow and transport in the natural environment: Advances and applications 183–229.
  * Flack & Schultz (2010) Flack, Karen A & Schultz, Michael P 2010 Review of hydraulic roughness scales in the fully rough regime. J. Fluids Eng. 132 (4).
  * Fourrière et al. (2010) Fourrière, A., Claudin, P. & Andreotti, B. 2010 Bedforms in a turbulent stream: formation of ripples by primary linear instability and of dunes by nonlinear pattern coarsening. J. Fluid Mech. 649, 287–328.
  * Frederick & Hanratty (1988) Frederick, K. A. & Hanratty, Thomas J. 1988 Velocity measurements for a turbulent nonseparated flow over solid waves. Exp. Fluids 6 (7), 477–486.
  * Gillette et al. (1998) Gillette, D.A., Marticorena, B. & Bergametti, G. 1998 Change in the aerodynamic roughness height by saltating grains: Experimental assessment, test of theory, and operational parameterization. J. Geophys. Res. Atmos. 103 (D6), 6203–6209.
  * Gillies et al. (2007) Gillies, J.A., Nickling, W.G. & King, J. 2007 Shear stress partitioning in large patches of roughness in the atmospheric inertial sublayer. Bound. Layer Meteorol. 122 (2), 367–396.
  * Gomé et al. (2022) Gomé, S., Tuckerman, L.S. & Barkley, D. 2022 Extreme events in transitional turbulence. Philos. Trans. Royal Soc. A 380 (2226), 20210036.
  * Gunn & Jerolmack (2022) Gunn, A. & Jerolmack, D.J. 2022 Conditions for aeolian transport in the solar system. Nat. Astron. 6 (8), 923–929.
  * Han et al. (2015) Han, C., Ma, Y., Su, Z., Chen, X., Zhang, L., Li, M. & Sun, F. 2015 Estimates of effective aerodynamic roughness length over mountainous areas of the tibetan plateau. Q. J. R. Meteorol. Soc. 141 (689), 1457–1465.
  * Hanratty (1981) Hanratty, T.J. 1981 Stability of surfaces that are dissolving or being formed by convective diffusion. Annu. Rev. Fluid Mech. 13 (1), 231–252.
  * Hobson et al. (1999) Hobson, J.M., Wood, N. & Brown, A.R. 1999 Large-eddy simulations of neutrally stratified flow over surfaces with spatially varying roughness length. Q. J. R. Meteorol. Soc. 125 (558), 1937–1958.
  * Hunt et al. (1988) Hunt, J.C.R., Leibovich, S. & Richards, K.J. 1988 Turbulent shear flows over low hills. Q. J. R. Meteorol. Soc. 114 (484), 1435–1470.
  * Jackson & Hunt (1975) Jackson, P.S. & Hunt, J.C.R. 1975 Turbulent wind flow over a low hill. Q. J. R. Meteorol. Soc. 101 (430), 929–955.
  * Jacobs (1989) Jacobs, S.J. 1989 Effective roughness length for turbulent flow over a wavy surface. J. Phys. Oceanogr. 19 (7), 998–1010.
  * Jasinski & Crago (1999) Jasinski, M.F. & Crago, R.D. 1999 Estimation of vegetation aerodynamic roughness of natural regions using frontal area density determined from satellite imagery. Agric. For. Meteorol. 94 (1), 65–77.
  * Jia et al. (2017) Jia, P., Andreotti, B. & Claudin, P. 2017 Giant ripples on comet 67P/Churyumov–Gerasimenko sculpted by sunset thermal wind. Proc. Natl. Acad. Sci. U.S.A. 114 (10), 2509–2514.
  * Jiménez (2004) Jiménez, J. 2004 Turbulent flows over rough walls. Annu. Rev. Fluid Mech. 36, 173–196.
  * Johnson & Johnston (1989) Johnson, P.L. & Johnston, J.P. 1989 Active and inactive motions in a turbulent boundary layer-interactions with free-stream turbulence. In 7th Symposium on Turbulent Shear Flows 2, 20–22.
  * Kadivar et al. (2021) Kadivar, M., Tormey, D. & McGranaghan, G. 2021 A review on turbulent flow over rough surfaces: Fundamentals and theories. Int. J. Thermofluids. 10, 100077.
  * Kametani et al. (2015) Kametani, Y., Fukagata, K., Örlü, R. & Schlatter, P. 2015 Effect of uniform blowing/suction in a turbulent boundary layer at moderate reynolds number. Int. J. Heat Fluid Flow 55, 132–142.
  * Kamphuis (1974) Kamphuis, J.W. 1974 Determination of sand roughness for fixed beds. J. Hydraul. Res. 12 (2), 193–203.
  * Kent et al. (2017) Kent, C.W., Grimmond, S. & Gatey, D. 2017 Aerodynamic roughness parameters in cities: Inclusion of vegetation. J. Wind Eng. Ind. Aerodyn. 169, 168–176.
  * King et al. (2006) King, J., Nickling, W.G. & Gillies, J.A. 2006 Aeolian shear stress ratio measurements within mesquite-dominated landscapes of the chihuahuan desert, new mexico, usa. Geomorphology 82 (3-4), 229–244.
  * Lancaster & Baas (1998) Lancaster, N. & Baas, A. 1998 Influence of vegetation cover on sand transport by wind: field studies at owens lake, california. Earth Surf. Process. Landf. 23 (1), 69–82.
  * Lapôtre et al. (2016) Lapôtre, M., Ewing, R.C., Lamb, M.P., Fischer, W.W., Grotzinger, J.P., Rubin, D.M., Lewis, K.W., Ballard, M.J., Day, M., Gupta, S. & Others 2016 Large wind ripples on mars: A record of atmospheric evolution. Science 353 (6294), 55–58.
  * Lapôtre et al. (2018) Lapôtre, M.G.A., Ewing, R.C., Weitz, C.M., Lewis, K.W., Lamb, M.P., Ehlmann, B.L. & Rubin, D.M. 2018 Morphologic diversity of martian ripples: Implications for large-ripple formation. Geophys. Res. Lett. 45 (19), 10–229.
  * Lee & Moser (2015) Lee, M. & Moser, R.D. 2015 Direct numerical simulation of turbulent channel flow up to ${Re}_{\tau}\approx 5200$. J. Fluid Mech. 774, 395–415.
  * Maaß & Schumann (1996) Maaß, C. & Schumann, U. 1996 Direct numerical simulation of separated turbulent flow over a wavy boundary. Flow Simulation with High-Performance Computers II: DFG Priority Research Programme Results 1993–1995.
  * Maat et al. (1991) Maat, N., Kraan, C. & Oost, W.A. 1991 The roughness of wind waves. Bound. Layer Meteorol. 54 (1), 89–103.
  * Marticorena & Bergametti (1995) Marticorena, B. & Bergametti, G. 1995 Modeling the atmospheric dust cycle: 1. design of a soil-derived dust emission scheme. J. Geophys. Res. Atmos. 100 (D8), 16415–16430.
  * Narteau et al. (2009) Narteau, C., Zhang, D., Rozier, O. & Claudin, P. 2009 Setting the length and time scales of a cellular automaton dune model from the analysis of superimposed bed forms. J. Geophys. Res. Earth. Surf. 114 (F3).
  * Nield et al. (2013) Nield, J.M., King, J., Wiggs, G.F.S., Leyland, J., Bryant, R.G., Chiverrell, R.C., Darby, S.E., Eckardt, F.D., Thomas, D.S.G., Vircavs, L.H. & Others 2013 Estimating aerodynamic roughness over complex surface terrain. J. Geophys. Res. Atmos. 118 (23), 12–948.
  * Nikuradse (1933) Nikuradse, J. 1933 Laws of flow in rough pipes. NACA Technical Memorandum 1292 .
  * Owen (1964) Owen, P.R. 1964 Saltation of uniform grains in air. J. Fluid Mech. 20 (2), 225–242.
  * Park & Choi (1999) Park, J. & Choi, H. 1999 Effects of uniform blowing or suction from a spanwise slot on a turbulent boundary layer flow. Phys. Fluids 11 (10), 3095–3105.
  * Paul-Limoges et al. (2013) Paul-Limoges, E., Christen, A., Coops, N.C., Black, T.A. & Trofymow, J.A. 2013 Estimation of aerodynamic roughness of a harvested douglas-fir forest using airborne lidar. Remote Sens. Environ. 136, 225–233.
  * Perrard et al. (2019) Perrard, S., Lozano-Durán, A., Rabaud, M., Benzaquen, M. & Moisy, F. 2019 Turbulent windprint on a liquid surface. J. Fluid Mech. 873, 1020–1054.
  * Pope (2000) Pope, S.B. 2000 Turbulent flows. Cambridge university press.
  * Prandtl (1925) Prandtl, L. 1925 Bericht über untersuchungen zur ausgebildeten turbulenz. Z. Angew. Math. Mech. 5 (2), 136–139.
  * Prigent et al. (2005) Prigent, C., Tegen, I., Aires, F., Marticorena, B. & Zribi, M. 2005 Estimation of the aerodynamic roughness length in arid and semi-arid regions over the globe with the ers scatterometer. J. Geophys. Res. Atmos. 110 (D9).
  * Raupach (1991) Raupach, M.R. 1991 Saltation layers, vegetation canopies and roughness lengths. In Aeolian Grain Transport 1, 83–96. Springer.
  * Raupach (1992) Raupach, M.R. 1992 Drag and drag partition on rough surfaces. Bound. Layer Meteorol. 60 (4), 375–395.
  * Raupach et al. (1991a) Raupach, M.R., Antonia, R.A. & Rajagopalan, S. 1991a Rough-wall turbulent boundary layers. Appl. Mech. Rev. 44 (1), 1–25.
  * Raupach et al. (1993) Raupach, M.R., Gillette, D.A. & Leys, J.F. 1993 The effect of roughness elements on wind erosion threshold. J. Geophys. Res. Atmos. 98 (D2), 3023–3029.
  * Raupach et al. (1991b) Raupach, M. R., Antonia, R. A. & Rajagopalan, S. 1991b Rough-Wall Turbulent Boundary Layers. Applied Mechanics Reviews 44 (1), 1–25.
  * Richards (1980) Richards, K.J. 1980 The formation of ripples and dunes on an erodible bed. J. Fluid Mech. 99 (3), 597–618.
  * Rubanenko et al. (2022) Rubanenko, L., Lapôtre, M.G.A., Ewing, R.C., Fenton, L.K. & Gunn, A. 2022 A distinct ripple-formation regime on mars revealed by the morphometrics of barchan dunes. Nat. Commun. 13 (1), 7156.
  * Sadique et al. (2017) Sadique, J., Yang, X.I.A., Meneveau, C. & Mittal, R. 2017 Aerodynamic properties of rough surfaces with high aspect-ratio roughness elements: effect of aspect ratio and arrangements. Bound. Layer Meteorol. (2), 203–224.
  * Schlichting (1937) Schlichting, H. 1937 Experimental investigation of the problem of surface roughness. NACA Technical Memorandum 823 .
  * Schlichting & Gersten (2000) Schlichting, H. & Gersten, K. 2000 Boundary Layer Theory. 8th edn. Springer.
  * Schmid & Bünzli (1995) Schmid, H.P. & Bünzli, B. 1995 The influence of surface texture on the effective roughness length. Q. J. R. Meteorol. Soc. 121 (521), 1–21.
  * Schultz & Flack (2009) Schultz, M.P. & Flack, K.A. 2009 Turbulent boundary layers on a systematically varied rough wall. Phys. Fluids 21 (1), 015104.
  * Sherman & Farrell (2008) Sherman, D.J. & Farrell, E.J. 2008 Aerodynamic roughness lengths over movable beds: Comparison of wind tunnel and field data. J. Geophys. Res. Earth Surf. 113 (F2).
  * Smeets & Van den Broeke (2008) Smeets, C.J.P.P. & Van den Broeke, M.R. 2008 Temporal and spatial variations of the aerodynamic roughness length in the ablation zone of the greenland ice sheet. Bound. Layer Meteorol. 128 (3), 315–338.
  * Sykes (1980) Sykes, R.I. 1980 An asymptotic theory of incompressible turbulent boundary layer flow over a small hump. J. Fluid Mech. 101 (3), 647–670.
  * Taylor et al. (1976) Taylor, P.A., Gent, P.R. & Keen, J.M. 1976 Some numerical solutions for turbulent boundary-layer flow above fixed, rough, wavy surfaces. Geophys. J. Int. 44 (1), 177–201.
  * Taylor et al. (1989) Taylor, P.A., Sykes, R.I. & Mason, P.J. 1989 On the parameterization of drag over small-scale topography in neutrally-stratified boundary-layer flow. Bound. Layer Meteorol. 48 (4), 409–422.
  * Tsai et al. (2010) Tsai, J.L., Tsuang, B.J., Lu, P.S., Chang, K.H., Yao, M.H. & Shen, Y. 2010 Measurements of aerodynamic roughness, bowen ratio, and atmospheric surface layer height by eddy covariance and tethersonde systems simultaneously over a heterogeneous rice paddy. J. Hydrometeorol. 11 (2), 452–466.
  * Tsuang et al. (2003) Tsuang, B.J., Tsai, J.L., Lin, M.D. & Chen, C.L. 2003 Determining aerodynamic roughness using tethersonde and heat flux measurements in an urban area over a complex terrain. Atmos. Environ. 37 (14), 1993–2003.
  * Tuckerman et al. (2020) Tuckerman, L.S., Chantry, M. & Barkley, D. 2020 Patterns in wall-bounded shear flows. Annu. Rev. Fluid Mech. 52, 343–367.
  * Valance et al. (2015) Valance, A., Rasmussen, K.R., El Moctar, A.O. & Dupont, P. 2015 The physics of aeolian sand transport. Comptes Rendus. Physique 16 (1), 105–117.
  * Van Driest (1956) Van Driest, E.R. 1956 On turbulent flow near a wall. J. Aeronaut. Sci. 23 (11), 1007–1011.
  * Van Rijn (1982) Van Rijn, L.C. 1982 Equivalent roughness of alluvial bed. J. Hydraul. Div. 108 (10), 1215–1218.
  * Van Rijn (1984) Van Rijn, L.C. 1984 Sediment transport, part iii: bed forms and alluvial roughness. J. Hydraul. Eng. 110 (12), 1733–1754.
  * Venditti et al. (2005) Venditti, J.G., Church, M. & Bennett, S.J. 2005 Morphodynamics of small-scale superimposed sand waves over migrating dune bed forms. Water Resour. Res. 41 (10), W10423.
  * Wei & Willmarth (1989) Wei, T. & Willmarth, W.W. 1989 Reynolds-number effects on the structure of a turbulent channel flow. J. Fluid Mech. 204, 57–95.
  * Wiberg & Nelson (1992) Wiberg, P.L. & Nelson, J.M. 1992 Unidirectional flow over asymmetric and symmetric ripples. J. Geophys. Res. Oceans 97 (C8), 12745–12761.
  * Wolfe & Nickling (1993) Wolfe, S.A. & Nickling, W.G. 1993 The protective role of sparse vegetation in wind erosion. Prog. Phys. Geogr. 17 (1), 50–68.
  * Xue et al. (2002) Xue, X., Wang, T., Sun, Q. & Zhang, W. 2002 Field and wind-tunnel studies of aerodynamic roughness length. Bound. Layer Meteorol. 104 (1), 151–163.
  * Zanke et al. (2022) Zanke, U., Roland, A. & Wurpts, A. 2022 Roughness effects of subaquaeous ripples and dunes. Water 14 (13), 2024.
  * Zilker et al. (1977) Zilker, D.P., Cook, G.W. & Hanratty, T.J. 1977 Influence of the amplitude of a solid wavy wall on a turbulent flow. part 1\. non-separated flows. J. Fluid Mech. 82 (1), 29–51.
  * Zilker & Hanratty (1979) Zilker, D.P. & Hanratty, T.J. 1979 Influence of the amplitude of a solid wavy wall on a turbulent flow. part 2. separated flows. J. Fluid Mech. 90 (2), 257–271.