# Linear models of strip-type roughness

D. Lasagna1<EMAIL_ADDRESS>G. Zampino2  B. Ganapathisubramani1
1Aeronautics and Astronautics, Faculty of Engineering and Physical Sciences,
University of Southampton, Hampshire, SO17 1BJ, UK 2 FLOW, Faculty of
Engineering Mechanics, KTH Royal Institute of Technology Stockholm, Sweden

###### Abstract

Prandtl’s secondary flows of the second kind generated by laterally-varying
roughness are studied in this work using the linearised Reynolds-Averaged
Navier-Stokes approach recently proposed in Zampino et al. (2022). The
momentum equations are coupled to the Spalart-Allmaras transport model for the
turbulent viscosity and the surface roughness is captured by adapting
established strategies for homogeneous rough walls to spanwise inhomogeneous
surfaces. Linearisation of the governing equations and of the roughness model
yields a linear framework that elucidates the linear mechanisms that play a
role in the generation of secondary flows. In addition, the framework allows a
rapid exploration of the parameter space associated with heterogeneous rough
surfaces. Channel flow is considered, with high and low roughness strips
arranged symmetrically on the two channel walls. The strip width is varied
systematically to ascertain the role played by the length scale characterising
the heterogeneity on the size and intensity of secondary flows. It is found
that linear mechanisms, i.e. whereby secondary flows are interpreted as the
output response of the turbulent mean flow subjected to a forcing localised at
the wall, may be sufficient to explain such a role. In fact the model predicts
that secondary flows are most intense when the strip width is about 0.7 times
the half-channel height, in excellent agreement with available data. Further,
a unified framework to analyse combinations of heterogeneous roughness
properties and laterally-varying topographies, common in engineering
applications, is discussed. This framework yields two separate secondary-flow
inducing source mechanisms, i.e. the lateral variation of the virtual origin
from which the turbulent structure is allowed to develop and the lateral
variation of the streamwise velocity slip, capturing the
acceleration/deceleration perceived by the bulk flow over troughs and crests
of non-planar topographies.

###### keywords:

Linearised RANS equations, secondary flows, strip-type roughness

##  1 Introduction

Prandtl’s secondary flows of the second kind (Prandtl, 1952) emerge when a
turbulent flow develops over an heterogeneous surface with a lateral variation
of its properties. Two equivalent standpoints, based on the analysis of the
Reynolds-averaged equations, explain the formation of such currents. One
standpoint considers secondary currents as the product of the imbalance
between production and dissipation of turbulent kinetic energy induced by the
roughness heterogeneity (Hinze, 1973), whereby turbulence-rich fluid is
advected towards low-turbulence regions. The second standpoint considers the
streamwise vorticity balance (Perkins, 1970), whereby cross-stream gradients
of the Reynolds stress arising from the cross-stream velocity components
induce a turbulent torque that acts as a source term in the streamwise
vorticity equation (Castro & Kim, 2024). Overall, such mechanisms produce
large-scale counter-rotating longitudinal rolls appearing in the time-averaged
wall-bounded flow. The associated upwelling and downwelling motions produced
by the rolls induce a lateral distortion of the boundary layer height (Barros
& Christensen, 2014), together with alternating high and low streamwise
momentum regions (Mejia-Alvarez & Christensen, 2013; Willingham et al., 2014;
Anderson et al., 2015), arranged analogously to the classical roll-streak
pattern in shear flows (Brandt, 2014).

Secondary flows are commonly observed in many industrial and environmental
applications, where surfaces are either characterised by lateral variations of
the topography, i.e. the elevation, or of the friction, e.g. by means of
varying roughness properties. These two types of heterogeneity have been
idealised in the literature as ridge-type and strip-type roughness
configurations, respectively. The first type consists of longitudinal ribs
located on the otherwise smooth, planar surface, and having rectangular or
more complex cross-sections (Goldstein & Tuan, 1998; Hwang & Lee, 2018;
Zampiron et al., 2020; Castro et al., 2021; Long et al., 2023; Zampino et al.,
2023; Zhdanov et al., 2024), or alternatively smooth sinusoidal modulations of
the wall (Wang & Cheng, 2006; Vidal et al., 2018). The second type, the focus
of this work, consists of alternating regions of high and low friction.
Secondary motions over such an arrangement have been extensively characterised
experimentally (Bai et al., 2018; Wangsawijaya et al., 2020; Wangsawijaya &
Hutchins, 2022; Frohnapfel et al., 2024), and in numerical simulations
(Willingham et al., 2014; Anderson et al., 2015; Chung et al., 2018; Stroh et
al., 2020; Forooghi et al., 2020; Neuhauser et al., 2022; Schäfer et al.,
2022).

Despite the burgeoning interest in these flows and the intense examination of
their characteristics, there is a number of aspects clearly documented in the
literature for which a physics-based, mechanistic model is not available. The
first aspect is related to the marked dependence of the size and intensity of
secondary flows on one or more spanwise length scales characterising the
surface heterogeneity (Wangsawijaya et al., 2020; Yang & Anderson, 2018;
Medjnoun et al., 2018). For strip-type roughness this length scale is usually
expressed by the width $S$ of the strips. Consensus has emerged on the
existence of three separate regimes as $S$ varies in relation with the average
boundary layer thickness $\delta$ (Chung et al., 2018). When the strip width
is much smaller than the boundary layer thickness, $S\ll\delta$, secondary
flows are confined to the vicinity of the surface and do not strongly
influence the outer region. Conversely, when $S\gg\delta$ secondary flows are
localised in regions where the surface properties vary more rapidly, and wide
areas of local flow homogeneity are observed away from such regions
(Wangsawijaya et al., 2020). When $S\approx\delta$, the secondary flows are
most intense and can significantly influence the flow structure. Nevertheless,
a model that captures the nature of these regimes and identifies boundaries
between them is not available at present. In addition, most studies have
considered strips of equal width, but the width ratio between high and low
roughness strips is certainly important, as it is for rectangular ridges where
the ratio of recessed and elevated area influences the flow structure
(Medjnoun et al., 2020; Zampino et al., 2022).

A second aspect is related to the occurrence of tertiary flows, i.e. weak roll
structures adjacent to the dominant rolls, often associated to a reversal of
the vertical flow direction at the centre of the high (or low) roughness
strip, or at the centre of the ridge (or trough) (Vanderwel &
Ganapathisubramani, 2015). Tertiary flows are commonly observed over surfaces
with longitudinal ridges (e.g. Medjnoun et al. (2020)), especially when the
width of the troughs or of the ridges is large enough to accommodate multiple
streamwise vortices next to each other. For heterogeneous rough surfaces,
however, tertiary flows have not been observed. In fact, for wide strips,
cross-stream motions have been observed to be mostly confined in a roughly
square region around the transition between the strips. This applies to both
boundary layer (Wangsawijaya et al., 2020) and channel flows (Chung et al.,
2018; Neuhauser et al., 2022). One explanation may be that tertiary flows over
roughness strips might be difficult to discern in the mean flows obtained from
experiments or simulations, especially when instantaneous structures meander
quite significantly in the longitudinal direction (Zampiron et al., 2021;
Kevin et al., 2019), smearing weak cross-stream motions. Neuhauser et al.
(2022) hypothesised that the boundary conditions utilised in numerical
simulations to capture the roughness effect may also play a role, although
this hypothesis does not appear to explain why tertiary flows are not seen in
experiments.

A third aspect that still lacks a robust mechanistic explanation is motivated
by features of realistic surfaces in engineering and natural applications,
whereby lateral changes of the roughness height are almost invariably
accompanied by a lateral change in the elevation (Stroh et al., 2020; Schäfer
et al., 2022). Decoupling these two effects may be easier in numerical
simulations where the roughness heterogeneity is modelled by suitable spanwise
heterogeneous boundary conditions applied to an otherwise planar boundary of
the numerical domain (Chung et al., 2018; Neuhauser et al., 2022), but
requires care when setting up experiments with, e.g. sandpaper strips or in
roughness-resolving numerical simulations (Frohnapfel et al., 2024). One
explicit attempt to study the coupling between these two effects was carried
out by Stroh et al. (2020) and then later by Schäfer et al. (2022) who
performed a series of direct numerical simulations over surfaces characterised
by alternating rough and smooth regions. In their paper, Stroh et al. (2020)
completely resolved the surface roughness using an immersed boundary method
and studied three different configurations: the mean roughness height is (i)
lower, (ii) equal to, and (iii) higher than the elevation of the smooth
surface. The authors observed a change in the flow organisation moving from
case (i) to (iii) and vice-versa. This behaviour was not reproduced in more
recent roughness-resolving simulations (Frohnapfel et al., 2024), which was
attributed to the importance of the strip width, relative to the roughness
height.

One last aspect for which a model does not seem to be available concerns the
relation between naturally-occurring Very-Large-Scale-Motions (VLSMs),
populating the log-layer over homogeneous surfaces, and secondary flows (Chung
et al., 2018; Lee et al., 2019; Wangsawijaya & Hutchins, 2022). It has been
speculated that secondary flows may be interpreted as VLSMs locked in place by
the surface heterogeneity, given some similarity in their features. This can
readily explain why secondary flows are most intense when $S\approx\delta$,
because the strip width is commensurate with the spanwise length scale of such
motions. Evidence shows that VLSMs and secondary flows do indeed coexist and
do interact to a significant extent, since energy from the former appear to
leak into the latter (Zampiron et al., 2020; Wangsawijaya et al., 2020).
However, the specific mechanism for which large scale structures residing in
the outer layer should be locked in place so effectively by the roughness
heterogeneity is not fully clear.

In recent work (Zampino et al., 2022), we developed a predictive framework to
understand how far can linear mechanisms go in explaining these aspects,
focusing on ridge-type roughness. The framework originates from the long line
of work that relies on the analysis of the Reynolds-Averaged Navier-Stokes
equations, augmented with a turbulent viscosity model and linearised about the
turbulent mean, to explain the structure of wall-bounded turbulence (see del
Álamo & Jiménez (2006); Pujals et al. (2009) and references therein).
Differently to previous efforts (Meyers et al., 2019), the framework utilises
the Spalart-Allmaras equation (Spalart & Allmaras, 1994) to capture turbulent
viscosity transport phenomena in combination with the nonlinear Quadratic
Constitutive Relations (QCR, Spalart (2000)) to model the anisotropy of the
Reynolds stress tensor, required to produce secondary motions (see Speziale
(1982); Speziale et al. (1991); Bottaro et al. (2006)). One first advantage of
the linear framework is its computationally efficiently. It thus enables the
vast parameter space characterising heterogeneous surfaces to be explored
rapidly, for instance to unravel the effect of ridge geometry (Zampino et al.,
2023). A second key advantage of the framework is that it provides a
perspective of secondary motions as being the output response of the turbulent
mean, subjected to a steady perturbation localised at the wall produced by the
surface heterogeneity. Linear input-output analysis of the Reynolds-Averaged
Navier-Stokes equations has demonstrated to be a powerful tool to explain the
structure of wall-bounded turbulence (Hwang & Cossu, 2010; McKeon & Sharma,
2010).

In this paper, we bring the framework introduced in Zampino et al. (2022) to
bear on the problem of strip-type roughness. The effect of the surface
roughness is introduced following well-established modelling strategies for
rough walls (Aupoix, 2007; Prakash & Laurendeau, 2020). In a nutshell, such
strategies consists of modifying the virtual origin from which the turbulent
flow develops, such as to obtain the desired shift of the logarithmic profile.
In practice, this requires using inhomogeneous boundary conditions for the
turbulent quantities and altering the distance from the wall for turbulence
models in which such a field variable is key, such as in the SA model. The
overall aim of this paper is to characterise the formation and structure of
secondary flows developing above strip-type roughness by means of the proposed
linear framework. This will allow fundamental insight into the linear
mechanisms that control such flows to be generated. Specifically, we apply the
proposed linear framework to flows in channels, and examine the role that the
surface arrangement plays on: (i) the strength of secondary motions as a
function of the strip width, identifying the three regimes discussed in (Chung
et al., 2018), (ii) the occurrence of tertiary flows as the duty cycle is
varied, (iii) the structure of low- and high-momentum pathways and finally
(iv) the combination of roughness and surface elevation effects.

The modelling framework, and its extension to rough surfaces, is presented in
section 2. Results are then reported in section 3. In section 4, the framework
is generalised to more complex surface heterogeneities, combining the effects
of roughness and surface elevation. Finally, conclusions are summarised in
section 5.

##  2 Methodology

###  2.1 Governing equations

The incompressible flow of a fluid with kinematic viscosity $\nu$ and density
$\rho$ is considered in a pressure driven channel with half-height $h$ and
subjected to a streamwise pressure gradient $\Pi$. The friction velocity
$u_{\tau}=\sqrt{\tau_{w}/\rho}$, with $\tau_{w}=h\Pi$ the mean friction,
yields the friction Reynolds number $\Rey_{\tau}=u_{\tau}h/\nu$. Index
notation is used for the Cartesian coordinates $x_{i}$ and velocities
components $u_{i}$. Lengths are generally normalised by $h$ and velocities by
$u_{\tau}$. The superscript $(\cdot)^{+}$ is omitted in the following to
reduce clutter, unless necessary.

Figure 1: Sketch of the bottom half of the channel with high and low roughness
strips, of width $S_{h}$ and $S_{l}$ respectively, placed symmetrically on the
bottom and upper (not shown) walls and aligned with the streamwise direction
$x_{1}$. This pattern repeats laterally with period $\Lambda=S_{h}+S_{l}$.
These dimensions are normalised with the channel half-height $h$. The
coordinate system is anchored at the bottom-plane and at the centre of one of
the high roughness strips. Owing to the symmetry of the strip arrangement on
the two walls, the shaded red area defines the fundamental repeating flow unit
in which the flow structure is visualised later in the paper.

The channel walls are covered by alternating strips of high and low roughness
having width $S_{h}$ and $S_{l}$, respectively, as shown in figure 1. The
strips are streamwise-aligned and are placed symmetrically on the two walls.
The pattern repeats with spanwise periodicity $\Lambda=S_{h}+S_{l}$, the
fundamental length scale. We also introduce the duty cycle $DC=S_{h}/\Lambda$
to characterise the relative width of the strips, and refer to $S$ as the
strip width when $S_{h}=S_{l}$, i.e. for $DC=0.5$.

The continuity and momentum equations are Reynolds averaged and made non-
dimensional using $h$, $u_{\tau}$ and $\rho$. Average and fluctuation
quantities are denoted by an overbar and a prime. For streamwise-aligned
strips, we assume a streamwise-independent time-averaged flow, i.e.
$\partial(\cdot)/\partial x_{1}\equiv 0$, which filters out the meandering of
secondary currents Zampiron et al. (2020). As a result, the mean pressure can
be eliminated by considering the mean streamwise vorticity equation and
introducing the streamfunction $\overline{\psi}$, satisfying
$\nabla^{2}\overline{\psi}=\overline{\omega}_{1}$ with

$\overline{\omega}_{1}={\frac{\partial\overline{u}_{3}}{\partial
x_{2}}-\frac{\partial\overline{u}_{2}}{\partial x_{3}}}$ (1)

the mean streamwise vorticity. The Reynolds-averaged equations for the
streamwise momentum and the streamfunction are then ∂¯ψ∂x2 ∂¯u1∂x3 \- & ∂¯ψ∂x3
∂¯u1∂x2 = 1+1τ ( ∂2¯u1∂x22 +  ∂2¯u1∂x32) + ∂τ12∂x2 +  ∂τ13∂x3,
∂2∂x2∂x3  [  (∂¯ψ∂x2)2 -( ∂¯ψ∂x3)2  ]  + ( ∂2∂x32 -  ∂2∂x22)∂¯ψ∂x2
∂¯ψ∂x3 =
1τ  ( ∂2∂x22 + ∂2∂x32 )2 ¯ψ + ∂2∂x2∂x3 ( τ33-τ22) + ( ∂2∂x22 -
∂2∂x32 ) τ23, where $\tau_{ij}=-\overline{u_{i}^{\prime}u_{j}^{\prime}}$ is
the Reynolds stress tensor.

###  2.2 Turbulence modelling

When the linear Boussinesq hypothesis is used to express the deviatoric
component of the Reynolds stresses as a function of the mean velocity
gradients, namely

$\tau_{ij}^{L}=2\nu_{t}S_{ij},$ (2)

where $\nu_{t}$ is the turbulent eddy viscosity and $S_{ij}$ is the mean
velocity gradient

$S_{ij}=\frac{1}{2}\left(\frac{\partial\overline{u}_{i}}{\partial
x_{j}}+\frac{\partial\overline{u}_{j}}{\partial x_{i}}\right),$ (3)

the Reynolds stresses in (1)b do not depend on the streamwise velocity. Then,
the streamfunction equation decouples from the streamwise momentum equation
and its solution is trivially $\overline{\psi}\equiv 0$, i.e. no secondary
flows are generated.

As extensively discussed in the literature (see e.g. Perkins (1970); Speziale
(1982); Bottaro et al. (2006)), a nonlinear stress model is needed to predict
Prandlt’s secondary flows of the second kind, produced by spatial gradients of
the anisotropy of the Reynolds stresses. Several approaches have been proposed
in literature (e.g. Speziale (1982); Chen et al. (1997)). Here we utilise the
Quadratic Constitutive Relation (QCR) nonlinear model presented in Spalart
(2000), whereby the deviatoric component of the Reynolds stresses becomes

$\tau_{ij}^{Q}=\tau_{ij}^{L}-c_{r1}\left[O_{ik}\tau_{jk}^{L}+O_{jk}\tau_{ik}^{L}\right],$
(4)

where $O_{ij}$ is the normalised rotation tensor defined as

$O_{ij}={2W_{ij}}/{\sqrt{\displaystyle{\frac{\partial\overline{u}_{m}}{\partial
x_{n}}\frac{\partial\overline{u}_{m}}{\partial x_{n}}}}},$ (5)

and $W_{ij}$ is the anti-symmetric part of the velocity gradient tensor. The
QCR model depends on a tuning single constant, whose value $c_{r1}=0.3$ was
calibrated to match the anisotropy of the outer region of wall-bounded
turbulent flows in Spalart (2000). The default value is used throughout the
paper.

To close the momentum equations, the Spalart-Allmaras (SA) turbulence model
(Spalart & Allmaras, 1994) is employed. The SA model is preferred here over
other commonly employed two-equation models because it can be linearised
relatively easily. In addition, the SA model was developed for attached shear
flows, hence it should provide satisfactory predictions for the present case.
The steady SA model defines a transport equation for the modified eddy
viscosity $\tilde{\nu}$, normalised with $u_{\tau}$ and $h$, which coincides
with the standard eddy viscosity everywhere except in the viscous sub-layer.
It is related to the turbulent viscosity by the relation
$\nu_{t}=\tilde{\nu}f_{v1}$, where $f_{v1}={\chi^{3}}/({\chi^{3}+c_{v1}^{3}})$
with $\chi=\Rey_{\tau}\tilde{\nu}$ and $c_{v1}$ a tuning constant. Its
transport equation is

$\overline{u}_{i}\frac{\partial\tilde{\nu}}{\partial
x_{i}}=c_{b1}\tilde{\mathcal{S}}\tilde{\nu}+\frac{1}{\sigma}\left\\{\frac{\partial}{\partial
x_{j}}\left[\left(\frac{1}{\Rey_{\tau}}+\tilde{\nu}\right)\frac{\partial\tilde{\nu}}{\partial
x_{j}}\right]+c_{b2}\frac{\partial\tilde{\nu}}{\partial
x_{j}}\frac{\partial\tilde{\nu}}{\partial
x_{j}}\right\\}-c_{w1}f_{w}\left(\frac{\tilde{\nu}}{d}\right)^{2},$ (6)

where the terms model, respectively, advection, production, diffusion and
destruction of the modified turbulent viscosity. The same model constants and
model definitions reported in Spalart & Allmaras (1994) are used in this work.

###  2.3 Roughness model for homogeneous surfaces

Many rough wall modelling strategies for homogeneous roughness (e.g., among
others, Durbin et al. (2000); Suga et al. (2006); Aupoix (2007); Brereton &
Yuan (2018); Prakash & Laurendeau (2020)) for RANS simulations rely on the
notion of equivalent sandgrain roughness $k_{s}^{+}$. These strategies,
described in this section, link the equivalent sandgrain roughness to suitable
non-zero turbulence quantities (the modified eddy viscosity for the SA model)
at the smooth, planar boundary of the numerical domain, to capture the
increased turbulence activity near the rough surface and obtain the desired
shift of the logarithmic profile as the main effect of the surface roughness.
No-slip boundary conditions are applied for the velocity. The mean turbulence
structure is assumed to develop from a new virtual origin, displaced beneath
the numerical boundary by a suitable distance $d_{0}^{+}$, to be determined
(Rotta, 1962). Relying on the outer-layer similarity hypothesis of Townsend
(1976), far away from the surface the law of the wall is preserved, and the
shift of the streamwise velocity profile observed over a rough surface is
captured by the empirical relation

$\left.\frac{\partial u_{1,r}}{\partial
x_{2}^{+}}\right|_{x_{2}^{+}}=\left.\frac{\partial u_{1,s}}{\partial
x_{2}^{+}}\right|_{x_{2}^{+}+d_{0}^{+}},$ (7)

where the subscripts $(\cdot)_{r}$ and $(\cdot)_{s}$ denote quantities over
the rough and smooth walls. Integrating this relation from the wall with
$u_{1,r}(x_{2}^{+}=0)=0$ yields

$u_{1,r}(x_{2}^{+})=u_{1,s}(x_{2}^{+}+d_{0}^{+})-u_{1,s}(d_{0}^{+}),$ (8)

which evaluated far away from the surface gives the logarithmic shift

$\Delta u_{1}=u_{1,r}(x_{2}^{+}\gg d_{0}^{+})-u_{1,s}(x_{2}^{+}\gg
d_{0}^{+})=u_{1,s}(d_{0}^{+})$ (9)

i.e. the distance $d_{0}^{+}$ can be found as the wall-normal coordinate where
the velocity over the smooth wall is equal to the desired velocity shift
$\Delta u_{1}$.

Then a roughness function that links the equivalent sandgrain roughness to the
shift of the log region is required. Among various options available, here we
use the Colebrook-Grigson roughness function (Grigson, 1992), given by

$\Delta
u_{1}=\frac{1}{\kappa}\log\left(1+\frac{k_{s}^{+}}{\exp(3.25\kappa)}\right).$
(10)

Although significant variations can be observed in the transitional regime,
little practical difference are evinced for the fully-rough regime in using
this and other models, such as Nikuradze’s roughness function (Aupoix, 2007).
Equations (10) and (9) allows the displacement $d_{0}^{+}$ to be expressed as
a function of the desired sandgrain roughness $k_{s}^{+}$.

With such information, a solution consistent with (7) can be found when the
eddy viscosity satisfies

$\nu_{t,r}(x_{2}^{+})=\nu_{t,s}(x_{2}^{+}+d_{0}^{+})$ (11)

This is achieved in two steps. First, the inhomogeneous boundary condition

$\tilde{\nu}_{r}(x_{2}^{+}=0)=\tilde{\nu}_{s}(d_{0}^{+})=d_{0}^{+}\kappa/\Rey_{\tau}$
(12)

is enforced to the modified eddy viscosity in the SA model, where the second
equality stems from the fact that $\tilde{\nu}$ varies linearly as $\kappa
x_{2}$ near the wall, where $\kappa$ is the von Kármán constant equal to 0.41.
Second, the distance $d$ between any point in the computational domain and the
nearest wall appearing in the SA model, as a fundamental field variable that
controls the balance between the production and destruction transport terms,
needs to be updated to reflect the location of the new virtual origin,
slightly below the numerical domain boundary. Overall, this procedure allows
the SA model to produce the desired shifted logarithmic velocity profile,
consistently with the updated boundary condition (12) on the modified eddy
viscosity.

Figure 2 illustrates example results at $\Rey_{\tau}=1000$, at which most of
the results presented in later sections were obtained.

Figure 2: Colebook-Grigson roughness function (a). Virtual origin $d_{0}^{+}$
as a function of the equivalent sandgrain roughness $k_{s}^{+}$ obtained by
using the smooth-wall SA velocity profile or the standard log-law (black
lines). The orange lines denote the derivative of $d_{0}^{+}$ with respect to
$k_{s}^{+}$ (axes on the right hand side). Numerical solutions of the
streamwise momentum equation and the SA turbulence model at order zero for
channels with smooth and homogeneous rough walls for $\Rey_{\tau}=1000$. Panel
(c) streamwise velocity, panel (d), turbulent viscosity. The circles in (a)
indicate the logarithmic shift associated to the profiles of panel (c),
calculated using the procedure in the text. For the logarithmic law, the
constants $\kappa=0.41$ and $A=5.1$ are used.

Solutions were obtained numerically with an in-house RANS code, based on a
Chebychev-collocation discretization method. The nonlinear system of algebraic
equations formed by the streamwise momentum equation and the SA equation was
solved using a Jacobian-free Newton–Krylov technique (Knoll & Keyes, ????),
using the “hookstep” approach of Viswanath (2007) to improve convergence.
Initial guesses for the streamwise velocity were obtained by first solving the
momentum equation using Cess’s analytical eddy viscosity profile (Reynolds &
Hussain, 1972). For a desired equivalent sandgrain roughness $k_{s}^{+}$ for
the homogeneous surface, the Colebrook-Grigson roughness function in panel (a)
of figure 2 is first used to obtain $\Delta u_{1}$. Using the smooth-wall
velocity profile obtained from the SA model, panel (c), the virtual origin
$d_{0}^{+}$ is obtained upon applying (9) (panel b). Repeating this procedure
for several equivalent sandgrain roughness yields the curve “SA-CG” in panel
(b). Clearly, the choice of the smooth-wall velocity profile influences the
results. For instance, coupling the Colebrook-Grigson formula the log-law
$u_{1,s}(x_{2}^{+})=\log(x_{2}^{+})/\kappa+A$, with $\kappa=0.41$ for
consistency with the standard SA model and $A=5.1$, yields

$d_{0}^{+}=\left(1+\frac{k_{s}^{+}}{\exp(3.25\kappa)}\right)\exp(-\kappa
A)\simeq 0.0326k_{s}^{+},$ (13)

denoted as “Log-CG” in the figure. This virtual origin is then used for the
boundary condition (12) and in the SA model. Overall, this yields shifted
velocity profiles, panel (c), that produce the desired $\Delta u_{1}$, denoted
by the circles in panel (a). It is worth pointing out that small variations of
the turbulent viscosity distribution in panel (d) are sufficient to produce
relatively significant alterations of the mean velocity profile. Numerically,
this makes the equations relatively stiff to solve.

###  2.4 Roughness model for heterogeneous surfaces

The equivalent sandgrain roughness is a dynamic parameter that is non
trivially related to the roughness geometry. For homogeneous roughness, it can
be readily estimated from correlations once the shift of the velocity profile
is known. However, for heterogeneous roughness, e.g. the present surface with
alternating strips, it is not immediately clear how one should assign an
equivalent sandgrain roughness to the two strips from velocity measurements,
as the flow structure and thus the logarithmic shift also depend on the
spatial distribution of the roughness properties (Wangsawijaya et al., 2020).
This conundrum is fundamentally the same discussed in numerical simulation
studies in which the roughness is not resolved but suitable boundary
conditions are applied at the smooth, planar boundary of the numerical domain.
In such strategies, the roughness heterogeneity can be modelled directly by a
lateral variation of the shear stress (Chung et al., 2018) or by a lateral
variation of the transversal slip length (Neuhauser et al., 2022). These
strategies require a model that links the boundary conditions to the desired
logarithmic shift. Such models are generally derived for homogeneous surfaces
and their applicability to heterogeneous surfaces may be questioned.

Here, given the lack of a better strategy, the just aforementioned approach is
adopted. Specifically, the alternating strips are defined by a spanwise
variation of the equivalent sandgrain roughness, following the expression

$k_{s}^{+}(x_{3})=k_{s}^{(0)}+\epsilon k_{s}^{(1)}(x_{3}),$ (14)

where $k_{s}^{(0)}$ is the reference, spatially-constant roughness height and
$k_{s}^{(1)}(x_{3})$ captures the variation of the roughness properties over
the two strips. For a unitary $h_{s}^{(1)}$ amplitude of the roughness
pattern, the spanwise variation is defined by the unitary peak-to-peak
amplitude, zero-mean function

$\displaystyle k_{s}^{(1)}(x_{3})=\begin{cases}1-DC&0\leq x_{3}\leq
S_{h}/2~{}\mathrm{and}~{}\Lambda-S_{h}/2\leq x_{3}\leq\Lambda\\\
-DC&S_{h}/2\leq x_{3}\leq\Lambda-S_{h}/2,\end{cases}$ (15)

where $DC=S_{h}/\Lambda$ is the duty cycle of the roughness modulation.
Because the difference in roughness between the two strips defined by
$k_{s}^{(1)}(x_{3})$ is unitary, the parameter $\epsilon$ in equation (14)
controls the actual difference in roughness between the two strips. This
definition, rather than specifying the roughness of the two strips, or
considering roughness strips separated by smooth regions (Wangsawijaya et al.,
2020), is motivated by the fact that we consider in the present analysis the
asymptotic limit when $\epsilon$ tends to zero, so that linearised equations
governing the response of the turbulent shear flow developing over a rough
surface to a small spanwise variation of the roughness properties can be
obtained.

###  2.5 Linearisation of the Reynolds-averaged equations and of the roughness
model

The streamwise momentum, streamfunction and SA equations form a coupled system
of three nonlinear partial differential equations that can be solved for any
desired strip configuration. However, when the difference between the
properties of the high and low roughness strips is small, in the limit when
$\epsilon\ll 1$, the resulting flow structure can be thought of as the sum of
the flow in a channel with the homogeneous reference roughness $k_{s}^{(0)}$
and a small perturbation, produced by the surface roughness heterogeneity
$k_{s}^{(1)}(x_{3})$ and capturing the heterogeneous flow structure of
secondary flows. This small perturbation obeys a set of linear equations which
is much easier to solve, and only captures linear input-output mechanisms. To
derive such equations, a generic time-averaged flow quantity
$\overline{q}(x_{2},x_{3})$ is first expanded as

$\overline{q}(x_{2},x_{3})\simeq
q^{(0)}(x_{2})+\epsilon{q}^{(1)}(x_{2},x_{3}).$ (16)

Substituting this expression for all mean quantities in the Reynolds-averaged
equations and in the SA equation, and taking terms at order zero in
$\epsilon$, leads to the nonlinear equations governing the flow over the
homogeneous rough surface. These are solved first, using the approach
discussed in section 2.3. At first order, the equations governing the
perturbation of the streamwise velocity and the streamfunction are -∂ψ(1)∂x3
Γ&=1τ(∂2∂x22+∂2∂x32) u1(1)+∂τ12(1)∂x2+∂τ13(1)∂x3,
0=1τ( ∂2∂x22+∂2∂x32)2 ψ(1)+∂2∂x2∂x3 ( τ33(1)-τ22(1))+(
∂2∂x22-∂2∂x32) τ23(1), where we define $\Gamma=\partial
u_{1}^{(0)}/\partial x_{2}$, showing that the zero-order solution, through the
mean velocity gradient $\Gamma$, needs to be available for the solution of the
first-order equations. The spanwise variation of the surface roughness and the
resulting secondary currents also introduce an alteration of the spatial
organisation of the turbulent viscosity through an alteration of the balance
of the transport terms in the SA equation. The linearised SA equation
governing such organisation is strongly coupled to the streamwise momentum and
the streamfunction equations and contributes to the perturbation of the
Reynolds stresses entering (2.5). A full analysis of the linearisation
procedure is reported in our previous work (Zampino et al., 2022), and is
omitted here for brevity.

It is worth noting that the streamfunction equation contains the perturbation
of the Reynolds stresses originating from the cross-stream velocity
components, as is well known (Perkins, 1970). Although at order zero, these
terms exhibit negligible influence, at order one the perturbation of the
Reynolds stress tensor becomes pivotal to couple the two equations in the
differential system (2.5). Here, the first-order stresses are found by
expanding the nonlinear Reynolds stress model (4) in a Taylor series in
$\epsilon$ and retaining the first order terms (see Zampino et al. (2022) for
definitions). After the linearisation, these stresses depend on both the
perturbation streamwise velocity and the perturbation streamfunction, leading
to a set of three coupled linear differential equations.

To obtain boundary conditions for the field variables, the wall roughness
treatment model discussed in section 2.4 needs to be linearised. The key idea
is that small spanwise perturbations of the equivalent sandgrain roughness are
modelled as small spanwise variations of the virtual origin. More formally,
over the heterogeneous surface given by (14), the shift of the virtual origin
varies according to

$d_{0}^{+}(x_{3})=d_{0}^{(0)}+\epsilon d_{0}^{(1)}(x_{3})$ (17)

where $d_{0}^{(0)}$ is the shift of the virtual origin of the reference
homogeneous surface with equivalent sandgrain roughness $k_{s}^{(0)}$. On the
other hand, the first-order term can be found by differentiating numerically
the curve reported in figure 2.3(b) at $k_{s}^{+}=k_{s}^{(0)}$ and using (15),
since

$d_{0}^{(1)}=\left.\frac{\mathrm{d}d_{0}^{+}}{\mathrm{d}k^{+}_{s}}\right|_{k_{s}^{(0)}}k_{s}^{(1)}.$
(18)

Asymptotically, for large $k_{s}^{(0)}$ and considering the log law, the
derivative of the curve in figure 2.3(b) tends to about 0.0326, implying that
a peak-to-peak variation of the equivalent sandgrain roughness of
$1/0.0326\approx 30$ is necessary to obtain a peak-to-peak variation of the
virtual origin equal to the viscous length scale. Once $d_{0}^{(1)}$ is known
linearising (12), the wall condition for the modified eddy viscosity

$\tilde{\nu}^{(1)}(x_{2}=0)=d_{0}^{(1)}\kappa/\Rey_{\tau},$ (19)

showing that the spanwise variation of the roughness properties is modelled as
a lateral change in the eddy viscosity at boundary of the numerical domain.
Finally, homogeneous boundary conditions are used for the streamwise velocity
perturbation and the streamfunction perturbation and its wall-normal
derivative.

One last remark is in order. A sensible question is whether the inclusion of
the linearised turbulence model to describe the perturbation of the turbulent
viscosity is really necessary, given that this is not customary in many
previous studies using linearised Navier-Stokes equations. On the one hand,
the strength of the mean flow response to a lateral perturbation of the
surface attributes is likely to depend in large part on the well-known
selective amplification properties of the linearised Navier-Stokes operator,
which are largest when such perturbation occurs at a specific spanwise length
scale, and much less on the linearised turbulence model. On the other hand,
the SA model provides a means to model realistically the effect of the surface
heterogeneity, because it provides clear insight into how the perturbation of
the effective distance $d$ influences the perturbation of the turbulent
viscosity field. As described in section 4, the lateral perturbation of the
distance $d$ is the dominant source mechanism that leads to secondary flows in
the present framework. In principle, one could first introduce an ansatz on
the perturbation of the turbulent viscosity and then only solve the linearised
Navier-Stokes equations. However, given the range of transport phenomena
modelled by the SA equation, defining the correct ansatz does not appear to be
a straightforward task.

###  2.6 Numerical solution of the linearised equation

The spanwise variation of the equivalent roughness height can be modeled as a
square wave approximated by the cosine series

$k_{s}^{(1)}=\sum_{n=1}^{\infty}\,k_{s}^{n}\cos\left(n\frac{2\pi}{\Lambda}x_{3}\right).$
(20)

The coefficients $k_{s}^{n}$ can be calculated analytically for each
combination of widths $S_{h}$ and $S_{l}$ and the corresponding coefficients
$d_{0}^{n}$ for the spanwise variation of the virtual origin $d_{0}^{(1)}$ are
found by using (18).

Expanding the unknown field variables at first-order in series, e.g. for the
streamwise velocity

$u_{1}^{(1)}(x_{2},x_{3})=\sum_{n=1}^{\infty}\,\hat{u}_{1}(x_{2};n)\cos\left(n\frac{2\pi}{\Lambda}x_{3}\right),$
(21)

and substituting these expressions in the linearised equations leads to one
set of three linear ordinary differential equations in $x_{2}$ for each
integer wavenumber. These are solved using a Chebyshev-collocation method.
Although the field variable $d$ in the SA model has a sharp cusp at $x_{2}=0$
and hence a spectral technique is not ideal for the solution of this problem,
we have observed that the numerical method is robust enough to provide
accurate results when a sufficiently fine grid is used. In the following
simulations, we used up to 252 collocation points. For the spanwise
discretisation, we observed that solutions converge relatively rapidly with
the number of Fourier modes retained in the expansion (20). This can be
motivated by the observation that, far away from the wall, only large-scale
perturbations of the surface features can influence the flow structure, while
the effect of small-scale perturbations, i.e. sharp gradients of the boundary
conditions (Neuhauser et al., 2022), decays more rapidly with the distance
from the wall (Meyers et al., 2019). In practice, we checked that results did
not change significantly when doubling the number of spanwise modes. The final
solution is then found by combining the solutions at each wavenumber, as the
superposition principle applies. One important implication of this property is
that the flow structure over surfaces with complex topographic/roughness
characteristics (Mejia-Alvarez & Christensen, 2013; Barros & Christensen,
2014) may be rationalised and better understood by decomposing the surface
forcing into its constitutive components. The strength of the mean flow
response at each spanwise length scale will then depend on the amplitude of
such components times a factor that captures the selective amplification
properties of the linearised Navier-Stokes operator (Chernyshenko & Baig,
2005), augmented by the linearised SA equation.

It is worth noting that the spatially-constant component at $n=0$ does not
appear in (20) because of the assumption that the spanwise variation of the
roughness height given by (15) is zero-mean. It does also not appear in the
solution (21), again because the linearity of the model implies that a
perturbation of the surface properties at wavenumber $n$ only produces a
distortion of the time averaged flow at the same wavenumber. A corollary of
this property is that the present model does not predict any change in mean
friction drag, the subject of several recent studies (Hutchins et al., 2023;
Frohnapfel et al., 2024). In fact, the spanwise-constant component
$\hat{u}_{1}(x_{2};0)$ and thus the perturbation of the bulk velocity computed
from this profile, is identically zero. The model does indeed capture the
spanwise modulation of the streamwise velocity distribution, i.e. high and low
momentum pathways (Barros & Christensen, 2014), but second order effects in
$\epsilon$, i.e. finite amplitude perturbations of the roughness height, are
necessary to capture the change in friction (Zampino, 2023).

##  3 Structure and strength of secondary currents

The volume averaged kinetic energy of the cross-sectional velocity components
$\mathcal{K}$, defined as

$\mathcal{K}=\frac{1}{4\Lambda}\\!\int_{0}^{2}\int_{0}^{\Lambda}\left[u_{2}^{(1)}(x_{2},x_{3})^{2}+u_{3}^{(1)}(x_{2},x_{3})^{2}\right]\,\mathrm{d}x_{3}\,\mathrm{d}x_{2},$
(22)

is used here to characterise the strength of the secondary flows. We also use
the streamfunction peak $\max_{x_{2},x_{3}}|\psi^{(1)}(x_{2},x_{3})|$ to
quantify the cross-stream flow rate associated with the vortices, as in other
studies (Vidal et al., 2018). Note that these variables are scaled with
$u_{\tau}$ and $h$. The solution of the linearised equations for a given strip
configuration can be obtained quite rapidly, which enables a rapid exploration
of the parameter space $(S_{h},S_{l})$. Results are reported in the two top
panels of figure 3 for $\Rey_{\tau}=1000$, and using $h_{s}^{(0)}=180$. Panels
(c-d) shows cuts along lines for three duty cycles as a function of the
fundamental length scale $\Lambda$.

Figure 3: Contours of the volume averaged kinetic energy of the cross-stream
velocities $\mathcal{K}$ (panel a) and the streamfunction peak value
$\max_{x_{2},x_{3}}|\psi^{(1)}|$ (panel b) as a function of the width of the
high- and low-roughness strips. The Reynolds number is $\Rey_{\tau}=1000$ and
$k_{s}^{(0)}=180$. Cases at constant spatial fundamental length scale
$\Lambda$ are identified by the dashed diagonal lines with negative slope.
Markers identify cases discussed later in the text. Panels (c-d) shows the
same two quantities for $DC=0.2,0.3$ and 0.5, as a function of $\Lambda$.

We note that, as in our previous application of these techniques to surfaces
with longitudinal ridges (Zampino et al., 2022), the results of the linearised
model become asymptotically Reynolds number independent for high Reynolds
numbers, somewhat supporting the weak Reynolds number dependence documented in
the literature (Wangsawijaya & Hutchins, 2022). Hence, the considerations
presented here can be, fundamentally, applied to higher Reynolds number flows
relevant to applications. Note also that flow variables, such as the velocity
or streamfunction perturbations, are computed in the present linear modelling
framework per unit variation of the equivalent sandgrain roughness height
$k_{s}^{(1)}$ (scaled in inner units) in analogy to what was described for
ridge-type roughness in Zampino et al. (2022) where the same quantities are
obtained per unitary ridge height (scaled in outer units). Given that
experiments on secondary flows over heterogeneous surfaces are often conducted
on roughness strips with a considerable difference in roughness properties
(Wangsawijaya et al., 2020; Chung et al., 2018), the numerical values reported
here will appear quite small. For graphical convenience, quantities are pre-
multiplied by a large factor, e.g. $10^{8}$ in 3(a), in the figures.

Regardless of the quantity used for measuring the strength of the secondary
currents, a peak is observed for $S_{h}=S_{l}\simeq 0.7$, corresponding to a
fundamental length scale $\Lambda\approx 1.4$, although the streamfunction
peaks slightly later. Further, the quantities in figure 3 are symmetric with
respect to the line $DC=0.5$, where the strength peaks. For $\Lambda\gtrsim
2.5$ two peaks are observed, located symmetrically with respect to the line
$DC=0.5$. Examination of the flow structure for some of these cases indicates
that the secondary flows observed over the high- and low-roughness regions for
a generic configuration $(S_{h},S_{l})$ are identical in strength but opposite
in flow direction when the width of the two strips is swapped. The strip width
at which secondary currents are most intense reflects previous observations.
For instance, Chung et al. (2018), using LES simulations, and Wangsawijaya et
al. (2020), using experiments above spanwise-alternating smooth and rough
strips, reported a maximum intensity for the secondary flows when the width of
the strips is comparable with the boundary layer thickness. In particular,
Wangsawijaya et al. (2020) found that the swirl strength was largest among the
cases they considered for $S_{h}=S_{l}=0.62$.

The contours of the kinetic energy around the peak region appear elongated
along the line $S_{h}+S_{l}=const$. As a result, the three cuts along lines at
constant duty cycle all display a peak for $\Lambda\simeq 1.4$, which may
interpreted as $\Lambda$ being the relevant length scale. This is partly
correct, as the response to sinusoidal perturbations of the sandgrain roughnes
does indeed peak for this length scale, as for ridge-type roughness (Zampino
et al., 2022). However, there is a marked effect of the relative size of the
strips away from the peak and the two width are indeed necessary to correctly
characterise the response. The maps of figure 3 show strong similarities with
the maps displayed in Zampino et al. (2022) and Zampino et al. (2023) for the
ridge-type roughness, as the peak amplification occurs for similar values of
$\Lambda$. This similarity suggests that the selective amplification of
secondary flows is an intrinsic property of the mean flow, and perhaps less
strongly an effect of the type of forcing, e.g. whether it is produced by
elevation (ridges) or roughness (strips) variations. In this regard, there has
been recent discussion on the relation and co-existence between secondary
currents and very large scale motions (VLSM, e.g. Lee et al. (2019)). One
speculation is that secondary currents are naturally-occurring VLSMs that are
phase locked spatially by the heterogeneous surface (Chung et al., 2018;
Wangsawijaya et al., 2020) and emerge in the time-averaged flow. The present
approach, based on the Reynolds-averaged equations where the concept of VLSMs
does not apply immediately, suggests that secondary currents may be
interpreted as the time-averaged response of a forcing localised near the wall
and produced by gradients of the turbulent stresses. The linearised Navier-
Stokes operator with its selective amplification properties then produces more
intense time-averaged structures at specific forcings wavelengths, when
$\Lambda\simeq 1.4$. This hypothesis leverages the same physics used in
transient growth analysis studies (del Álamo & Jiménez, 2006; Cossu et al.,
2009; Pujals et al., 2009) to explain the formation of coherent structures in
shear flows from properties of the Orr-Sommerfeld-Squire equations, except
that here we consider the steady response to a given steady perturbation
localised at the wall rather than the transient amplification of an optimal
initial perturbation.

###  3.1 Flow structure

To visualise the cross-stream structure of the secondary currents, fields of
the wall-normal and spanwise velocity perturbation are reported in figure 4,
for duty cycle $DC=0.5$, i.e. for equal length of the high and low roughness
strips, for a range of strip widths $S$.

Figure 4: Colour maps of the wall-normal velocity perturbation $u_{2}^{(1)}$
over roughness strips at $DC=0.5$ for several strip widths $S$ (left column).
The field in the fundamental domain (see figure 1) is repeated four and two
times for the two narrowest strip cases. Contour lines of the perturbation of
the streamfunction $\psi^{(1)}$ are also reported to better describe the
secondary flows. Negative $\psi^{(1)}$ are indicated by using dashed lines.
Contours of the spanwise velocity perturbation $u_{3}^{(1)}$ (right column).
The Reynolds number is $\Rey_{\tau}=1000$, for $h_{s}^{(0)}=180$. The
darker/lighter patches denote the high/low roughness strips.

These configurations corresponds to the star markers in figure 3(a). Contour
lines of the streamfunction are also reported. The roughness strips produce
two counter-rotating vortices inducing a down-welling over the high-roughness
regions and an up-welling over the low-roughness regions. For narrow strips,
even narrower than what shown here, the vortices are confined in the near wall
region and the flow appears homogeneous at distances from the wall larger than
the strip width. Increasing the strip width, the vortices grow with $S$ until
they occupy the full half-height of the channel. For $S\simeq 0.7$, the cross-
stream components and in particular the wall-normal component is more intense
than all other cases considered. This configuration identifies the peak
amplification region of the maps of figure 3. Increasing $S$ further, the
strength of the wall-normal motions decreases slightly, but the volume
averaged strength of such motions decreases much further as the flow structure
converges to an idealised wide-strip asymptotic limit where moderately intense
cross-stream motions are only found in the immediate vicinity of the
transition between strips, with fluid is at rest in the “homogeneous” regions
above the centres of the roughness strips. This explains the trend of the two
quantities in figure 3(c). While $\mathcal{K}$ is a volume averaged quantity
and decreases with $S$ for $S\gg 1$, the streamfunction peak is a local
quantity that measures the strength of the individual vortex cores. This
quantity shows a first peak for $S\simeq 0.75$, and eventually tends to an
asymptotic value (with a larger amplitude) for large $S$, characterising the
strength of the “isolated vortex” regime.

Near the transition between strips, the spanwise component is particularly
intense in the near-wall region (right panels in figure 4). The spanwise
velocity peak is more intense than the vertical velocity peak in agreement
with observations Frohnapfel et al. (2024), and is localised in the near wall
region. Given that the spanwise velocity obeys no-slip condition, this results
in very high streamwise vorticity localised at the transition between strips.
Analyses not reported here show that the wall-normal location of the spanwise
velocity peak scales in inner units when the Reynolds number is increased,
while its magnitude becomes $\Rey_{\tau}$ independent.

From a qualitative viewpoint, the flow structure predicted by the present
linearised model resembles previous experimental observations (Wangsawijaya et
al., 2020) and numerical simulations (Chung et al., 2018; Neuhauser et al.,
2022). One aspect of discussion concerns the spanwise location of the
streamwise-aligned vortices with respect to the alternating pattern of
roughness. In the present case, the vortices are symmetrically located above
the interface between the strips. This is because the linearity of the
governing equations preserves the symmetry that exists across the jump. By
contrast, the experiments of Wangsawijaya et al. (2020) show that the centres
of the vortices are typically found over the low-roughness region. The same
was predicted by the simulations of Chung et al. (2018), using an
inhomogeneous shear-stress boundary condition to model the roughness. The
puzzling aspect is that one would initially attribute the displacement of the
vortices to nonlinear convective effects not captured by the linear model, as
if the vortices were transported towards the low-roughness region by the
relatively intense spanwise velocities near the wall associated to the
streamwise vorticity field. However, such a displacement is not observed in
the simulations of Neuhauser et al. (2022) who applied a Navier slip boundary
condition for the spanwise velocity to model the roughness, or occurs in the
opposite direction in more-recent roughness-resolving simulations of submerged
roughness strips (Frohnapfel et al., 2024).

From a quantitative viewpoint, a comparison with published results is slightly
less straightforward given the particular setup considered in this paper. For
this purpose, we use the channel-flow simulations of Chung et al. (2018), at
$\Rey_{\tau}=590$. The roughness strips are modelled by setting the shear
stress to 50% more and 50% less than the average shear stress. Assuming that
the low and high roughness strips correspond to the smooth and rough wall
regions, respectively, Chung et al. (2018) estimates that the equivalent
sandgrain roughness of the roughness patches is $k_{s}^{+}=205$. With such
settings, they observe maximum wall-normal velocities that peak at $0.3\div
0.4u_{\tau}$ (see their figure 9a. To match the roughness properties of these
simulations, one needs to recall that the solution produced by the present
linear model is defined per unit variation of the equivalent sandgrain
roughness between the strips. From the results of figure 4, maximum velocities
on the order of $17\times 10^{-4}u_{\tau}$ are obtained, for the optimal width
$S$. Multiplying this value by $k_{s}^{+}=205$ we obtain velocities on the
order of $0.35u_{\tau}$, in very good quantitative agreement with the
numerical simulations. Given that the intensity of the cross-stream velocity
components characterises somehow the equilibrium between source/sink
mechanisms of the streamwise vorticity balance (Stroh et al., 2016; Castro &
Kim, 2024), the favourable agreement with simulations suggests that such
mechanisms are correctly captured by the linearised RANS model. However, the
maximum wall-normal velocity in Chung et al. (2018) is obtained for a
relatively wide strip, $S=1.57$, while the present model indicates that the
peak occurs at $S\approx 0.7$, and lower velocities are observed for $S=1.57$.
It is argued that this is not a Reynolds number effect, but it is due to the
vortices in Chung et al. (2018) being, as discussed, closer to each other than
the strip width $S$ would suggest, resulting in larger induced velocities.
Evidence for this is given by the fact that the wall-normal velocities depend
on the duty cycle $DC$, as shown later in figure 8, and peak when the vortices
are artificially pushed together by a narrow low-roughness strip,
$S_{L}<S_{h}$.

###  3.2 High- and low-momentum pathways

A unique characteristic of flows over heterogeneous surfaces is that the
longitudinal secondary currents are flanked by high- and low-momentum regions,
produced by the vertical “pumping” of high and low momentum fluid,
respectively, induced by the vortical motions (Mejia-Alvarez & Christensen,
2013; Willingham et al., 2014; Barros & Christensen, 2014), which determine a
significant spanwise alteration of the boundary layer depth. Although a wall-
bounded flow above an homogeneous rough surface displays a positive deficit of
the streamwise velocity due to the flow deceleration induced by the surface
roughness (Jiménez, 2004), the co-location between the high and low momentum
regions and the roughness strips is somehow counter-intuitive as faster flow
can be found for certain conditions above the high-roughness strips. In this
section, we demonstrate that the present framework clearly captures this
phenomenon.

Figure 5: Contours of the streamwise velocity perturbation $u_{1}^{(1)}$ over
roughness strips at $DC=0.5$ for several strip widths $S$, without QCR model
(left panels) and with QCR model (right panels). The field in the fundamental
domain (see figure 1) is repeated four and two times for the two narrowest
strip cases. Contour lines of the perturbation of the streamfunction
$\psi^{(1)}$ are also reported visualise the secondary flows. Negative
$\psi^{(1)}$ contours are indicated by dashed lines. The Reynolds number is
$\Rey_{\tau}=1000$, for $h_{s}^{(0)}=180$. The darker/lighter patches denote
the location of the high/low roughness strips.

In figure 5, contours of the streamwise velocity perturbation $u_{1}^{(1)}$
are shown for the same configurations of figure 4. This visualisation differs
to what customarily reported in previous work in that it shows the _velocity
deviation_ from the streamwise velocity distribution $u_{1}^{(0)}$ observed
over the homogeneous surface. The difference between the maps in the two
columns is that in the left panels secondary currents were artificially
“turned off” by setting the constant $c_{r1}$ in the nonlinear stress model to
zero, so that spanwise and wall-normal gradients of the mean streamwise
velocity do not produce any of the Reynolds stresses in the streamwise
vorticity equation necessary to sustain longitudinal vortices. In the right
panels, solutions for the standard value $c_{r1}=0.3$ are reported. Without
secondary currents, the flow experiences a net deceleration above the high-
roughness strips, especially in the near wall region. However, further away
from the wall, e.g. at the centre of the channel, the change in streamwise
velocity depends strongly on the strip width, because this parameter controls
the depth at which the roughness-induced deceleration “diffuses” in the shear
flow from the wall due to the turbulent viscosity field. For narrow strips,
the velocity deficit at the channel centre is small, and only becomes
significant when the the strip width is at least the half the half-height of
the channel. When secondary currents are “turned back on”, right panels of
figure 5, the streamwise velocity over the high-roughness strips is now
generally positive because of the down-welling motion in this region, and
vice-versa on the low-roughness strips. This only applies for $x_{2}\gtrsim
0.1$, because nearer to the wall the local roughness properties control
whether the flow is faster/slower. These motions produce dispersive stresses,
e.g. $u_{1}^{(0)}u_{2}^{(0)}$, that alter the equilibrium in the streamwise
direction and result in a non-trivial dependence of the streamwise velocity
from the wall. This, fundamentally, implies that the logarithmic velocity
distribution is significantly altered by the addition of the dispersive
stresses. However, as the strips become wider, secondary currents are not
intense enough to produce any significant alteration of the streamwise
momentum equilibrium and the deceleration effect produced by the high
roughness begins to dominate, starting from the region closest to the wall.
Overall, this is the same behaviour observed experimentally (Wangsawijaya et
al., 2020; Frohnapfel et al., 2024), where a downward/upward bulging of the
contours of the streamwise velocity are observed at the edge of the boundary
layer / near the wall.

The influence of the strip width on the perturbation of the streamwise
velocity field is summarised in figure 6. The streamwise velocity profile at
the centre of the high-roughness strip is extracted from several calculations
with $S$ in the range $[0.1,20]$, with and without the nonlinear Reynolds
stress model. These profiles are concatenated together to form the maps in
panels (a) and (b), respectively. For the case with $c_{r1}=0.3$, we also
report a similar plot for the wall-normal component, in panel (c). Given that
a logarithmic velocity shift is not a meaningful quantity to compute, we
report in panel (d) the velocity deviation in the centre of the channel, above
the centre of the high-roughness strips. In panel (c), the solid red line
denotes the wall-normal location of the perturbation streamfunction peak.

Figure 6: Streamwise velocity extracted on a vertical line at the centre of
the high-roughness strip for increasing strip width $S$, for $c_{r1}=0$, (a),
and the standard value 0.3, (b). Wall-normal velocity extracted on the same
vertical line for $c_{r1}=0.3$, (c). The dashed lines indicates wall-normal
distances growing linearly with $S$ for $DC=0.5$. Streamwise velocity
perturbation above the centre of the high-roughness strip, at the mid-plane,
for the same two cases as a function of the strip width, (d). The Reynolds
number is $\Rey_{\tau}=1000$ and $h_{s}^{(0)}=180$.

The three regimes discussed in Chung et al. (2018) can be clearly identified.
For wide strips, the velocity deficit tends to a value controlled by the
roughness function at $k_{s}^{(0)}$ for both cases, as the turbulent structure
over each strip tends to its equivalent over an homogeneous surface given that
the influence of neighbouring patches and of the secondary currents localised
at the transition between strips vanishes. In this regime, the model predicts
that the centre of the rolls is located at a distance from the wall of about
$0.42$, i.e. the rolls are space-filling in the vertical direction. The
streamwise velocity perturbation is also roughly constant as a function of the
distance from the wall for $S\gtrsim 8$, producing the expected shift of the
logarithmic velocity profile over each strip, regardless of whether the
nonlinear Reynolds stress is active or not.

For narrow strips, the height of the channel is much larger than the “depth”
at which the effect of the wall inhomogeneity is perceived, analogously to the
blending height concept for spatially varying roughness discussed in (Bou-Zeid
et al., 2007). For the case $c_{r1}=0$, without rolls, this depth (estimated
from the contours of the streamwise velocity perturbation) varies linearly
with $S$, as one would expect to see in a diffusion-driven problem, given that
the turbulent viscosity does not vary significantly except for near the wall.
For the case $c_{r1}=0.3$, with rolls, the depth measured in terms of the
wall-normal velocity also appears to increase linearly with $S$, at least for
the smallest $S$ analysed here. However, the linear scaling (the dashed lines)
appears to loose accuracy relatively rapidly, at $S\approx 0.25$, as soon as
the rolls occupy about half of the half-channel height. In between these two
regimes, in the “transitional regime” of Yang & Anderson (2018), the
streamwise velocity perturbation display a complex behaviour, highlighting the
significant lack of flow homogeneity. At $S=0.7$ the rolls are most intense
and induce the maximum perturbation at a wall-normal location close to the
centres of the rolls. However, the velocity perturbation reaches its peak at
the mid plane only at $S=0.95$ and then changes sign from $S\gtrsim 1.75$.
Interestingly, the model never predicts flow reversal when the strip width is
increased, in agreement with observations in the literature.

Overall, the present framework appears to capture correctly the three flow
regimes documented in the literature. The implication is that linear
mechanisms, whereby secondary flows may be interpreted as the output response
of the mean shear flow to a steady forcing localised at the wall, may be
sufficient to predict the size and strength of the rolls. Based on the
streamwise velocity evaluated at the mid-plane, the boundaries may be located
at $S\approx 0.4$ and $S\approx 8$, but differences may arise with alternative
criteria that consider the bulk of the flow.

###  3.3 Eddy viscosity perturbation

In many studies that have examined the properties of the Navier-Stokes
equation linearised around the turbulent mean profile and augmented with an
eddy viscosity to model the effect of the turbulent Reynolds stresses,
velocity perturbations resulting from optimal or stochastic forcing are
computed by assuming that the turbulent viscosity distribution is not affected
by the flow perturbation, and it is assumed a priori following analytical or
empirical distributions (Reynolds & Hussain, 1972; del Álamo & Jiménez, 2006;
Hwang & Cossu, 2010; Morra et al., 2019; Pickering et al., 2021). In the
present case, and unlike previous work, the governing equations must include a
transport equation for the turbulent viscosity. It is likely that the peak
location predicted by the present linear model (see figure 3) is purely
determined by the non-normality of the linearised Navier-Stokes operator, with
minimal effect of the turbulence model adopted. However, the inclusion of such
a model and the resulting perturbation of the turbulent viscosity
$\nu_{t}^{(1)}$ is key to capture the influence of the heterogeneous surface
roughness and the generation of secondary flows, as described in section 2.2.
In practice, spanwise gradients of the eddy viscosity produce, through the QCR
model, Reynolds stresses that act as source terms for the streamwise vorticity
equation, resulting in secondary motions.

Figure 7 shows contours of the perturbation turbulent viscosity for several
strip widths $S$. In the left panels, the QCR constant $c_{r1}$ has been set
to zero, to characterise the eddy viscosity distribution in the absence of
cross-stream motions and highlight the interaction with the streamwise
velocity fields of figure 5. In the right panels, the QCR constant is
$c_{r1}=0.3$. In general, positive eddy viscosity perturbations are observed
above the high roughness strips, and vice versa, reflecting the boundary
condition (19) and the altered distance $d$ from the wall. For narrow strips,
the eddy viscosity perturbation is confined near the wall, given that rapid
spatial variations of the eddy viscosity tend to be damped by the diffusion
term in the linearised SA equation. As the strip gets wider, more intense eddy
viscosity distributions are observed, reflecting the increased
acceleration/deceleration of the flow over the low and high roughness strips,
respectively.

Figure 7: Contours of the eddy viscosity perturbation $\nu_{t}^{(1)}$ over
roughness strips at $DC=0.5$ for several strip widths $S$, without QCR model
(left panels) and with QCR model (right panels). The field in the fundamental
domain (see figure 1) is repeated four and two times for the two narrowest
strip cases. Contour lines of the perturbation of the streamfunction
$\psi^{(1)}$ are also reported visualise the secondary flows. Negative
$\psi^{(1)}$ contours are indicated by dashed lines. The Reynolds number is
$\Rey_{\tau}=1000$, for $h_{s}^{(0)}=180$. The darker/lighter patches denote
the location of the high/low roughness strips.

Small, but likely significant, changes are observed when secondary flows are
“turned on”, right panels. It can be observed that the cross-stream motions
produce a further distortion of the eddy viscosity distributions. This is the
result of two mechanisms: a) the advection of turbulent viscosity of the
background flow operated by the vertical and lateral velocities and b) the
altered production of turbulent viscosity due to the alteration of wall-normal
and spanwise gradients of the streamwise velocity.

###  3.4 Tertiary structures and role of the duty cycle

Figure 8: Contours of the wall-normal (a-c) and streamwise (d-f) velocity
perturbation for increasing duty cycles, indicated in the figures, and
$\Lambda=2.8$. Panels (g) and (h) show the wall-normal and streamwise
velocity, respectively, extracted on the vertical line (shown in panels (a-f))
at the centre of the low-roughness strip $x_{3}/\Lambda=0.5$ for duty cycles
ranging from 0 (narrow high-roughness strip) to 1 (wide high-roughness
strips). The lines labeled (a-f) in panels (g, h) refer to the corresponding
panels. The Reynolds number is $\Rey_{\tau}=1000$ and $h_{s}^{(0)}=180$.

The present model shows that the duty-cycle appears to play an important role
in controlling the formation of tertiary structures, which have so far not
been observed for strip-type roughness (Neuhauser et al., 2022). This is shown
in figure 8, where contours of the wall-normal and streamwise perturbation
velocities are shown for increasing duty cycles, at $\Lambda=2.8$. These cases
correspond to the triangular marker in figure 3. For each of these solutions,
as well as for other solutions in the interval $DC=[0,1]$, the wall-normal and
streamwise velocity profiles are extracted at the centre of the low-roughness
strip. All these profiles are then combined together to form the colour map
reported in the bottom panels of figure 8. This visualisation suggests that
flow reversal over the low-roughness strip may begin to appear in practice as
the duty cycle is decreased to about $0.4$ and would be the most intense for
$DC\approx 0.23$ (i.e. for $S_{h}\approx 0.3S_{l}$). Incidentally, the
occurrence of flow reversal and tertiary structures explains the two “ears” of
the contours of the kinetic energy density map of figure 3(a), where the
cross-stream motions are relatively intense, compared to other duty cycles.
This analysis also suggests that the vertical velocities over the low-
roughness strip are most intense for a duty cycle equal to about 0.8 and not
for the symmetric case $DC=0.5$ which the kinetic energy density maps would
suggest. Arguably, this can be attributed to the constructive interference of
the wall-normal velocities induced by two neighbouring vortices, pushed closer
to each other by the decreasing width of the low-roughness strip. The
streamwise velocity field is also particularly affected by the duty cycle, as
faster or slower flow over the low roughness strip can be found depending on
the $DC$. The transport of fast/slow fluid operated by the cross-stream
velocities is particularly visible. For instance, for $DC=0.23$ the left
vortex, rotating counter-clockwise transports low-momentum fluid from the near
wall region at $x_{3}=0$ to its right flank, producing a negative velocity
streak at $x_{3}\approx 0.2$, and similarly for the other longitudinal vortex.

##  4 Generalising the framework to complex surface heterogeneities

It has been demonstrated that when conducting experiments (Wangsawijaya et
al., 2020) or roughness-resolving simulations (Stroh et al., 2020) over
realistic heterogeneous rough surfaces it is pivotal to ensure that the shear-
increasing effects of roughness are decoupled from the inevitable variation of
the mean surface height. Both roughness and elevation heterogeneity produce
secondary currents and therefore the combination of such effects can
significantly influence the strength and potentially the direction of the
resulting secondary motions (Schäfer et al., 2022; Frohnapfel et al., 2024).
In this section, we analyse this aspect through the lens of the linearised
model, to initiate the formulation of a unifying framework for flows over
complex heterogeneous surfaces.

###  4.1 Secondary-flow-inducing source mechanisms

In the present linearised framework, in the limit case where the spanwise
variation of the roughness or the elevation is small, ridge-type and strip-
type roughness are modelled with the same approach, whereby the flow-surface
interaction is captured through three separate source mechanisms corresponding
to three different inhomogeneous terms acting as forcing in the linearised
equations. To illustrate such source mechanisms, it is instructive to examine
ridge-type roughness considered in our previous work (Zampino et al., 2022)
where these mechanism all act simultaneously. In such case, the lateral
variation of the elevation was defined by a unitary peak-to-peak, zero-mean
function $f(x_{3})$ so that, e.g., the bottom wall of the channel is located
at $x_{2}=\epsilon f(x_{3})$ and the small parameter $\epsilon$ controls the
actual amplitude of the topography. The first source mechanism, denoted as $A$
in what follows, is mediated by the linearised boundary condition on the
streamwise velocity, e.g. on the lower wall,

$u_{1}^{(1)}(x_{2}=0)=-f(x_{3})\left.\frac{\partial u^{(0)}}{\partial
x_{2}}\right|_{x_{2}=0}=-f(x_{3})\Rey_{\tau}.$ (23)

Physically, this condition produces a velocity slip that captures the
acceleration and deceleration perceived by the bulk flow above the troughs and
the crests of the non-planar topography, respectively. What leads to the
formation of secondary flows is the resulting spanwise gradient of the
streamwise velocity in the near wall region. This gradient induce, via the QCR
model, spanwise gradients of anisotropic Reynolds stresses that then create
secondary flows. This mechanism does not seem to have been discussed
previously in the literature of secondary flows. More generally, we are not
aware of studies that consider the somewhat artificial surface arrangement
consisting of longitudinal, flush-mounted belts moving upstream/downstream
alternated with regions of solid wall. This setup would capture source
mechanism $A$ directly, and we predict that it could generate relatively
intense secondary currents with a down-welling over the upstream-moving belts.
When modelling strip-type roughness, this source mechanisms is not active
since no-slip boundary conditions for the streamwise velocity are used. This
is equivalent to stating that the mean height of the two roughness strips is
the same and the boundary of the numerical domain is at some suitable location
where the strip-averaged streamwise velocity goes to zero.

The second source mechanism, denoted as $B$, is mediated by the destruction
term of the SA transport equation, where the inverse of the squared distance
$d$ between a point in the numerical domain and the nearest “wall” models the
blocking effect of the wall (Spalart & Allmaras, 1994; Aupoix & Spalart,
2003). With such term, the SA model predicts an accurate log-layer, and thus
lateral perturbations of the distance $d$ induced by the elevation (via the
function $f(x_{3})$) or by the displacement of the virtual origin (via the
term $d_{0}^{(1)}(x_{3})$), produce a perturbations of the the log-layer and
spanwise gradients of the turbulent viscosity field. Crucially, the sign of
this source mechanism is opposite for strip-type and ridge-type roughness:
while locally increasing elevations correspond to locally reducing distances
$d$, increasing roughness produces increasing distances, as the virtual origin
is further displaced downwards beneath the boundary of the numerical domain.
In this regard, the framework suggests that the mean roughness height is not
the important factor. Rather, the displacement of the virtual origin, a
dynamic parameter that ultimately depends on the drag of the surface, is what
controls the intensity of the forcing and of the resulting secondary flows. As
a side note, while source mechanism $A$ acts as a boundary condition, source
mechanism $B$ acts at all wall distances, as it captures the perturbed
development of the wall-bounded flow from a different origin.

The third source mechanism, denoted as $B^{\prime}$, is mediated by the
inhomogeneous boundary condition on the modified turbulent viscosity. For
ridge-type roughness, this is

$\tilde{\nu}^{(1)}(x_{2}=0)=-f(x_{3})\left.\frac{\partial\tilde{\nu}^{(0)}}{\partial
x_{2}}\right|_{x_{2}=0}=-f(x_{3})\kappa.$ (24)

For strip-type roughness, the lateral variation of the virtual origin
$d_{0}^{(1)}(x_{3})$ plays the same role of the function $f(x_{3})$ describing
the ridge topography, but with an opposite effect as dictated by the boundary
condition (19). This condition must be applied consistently with source
mechanism $B$, so that the shifted eddy viscosity profile produced by such
mechanism is consistent with the boundary condition (24). Studies on modelling
roughness in RANS simulations (Aupoix & Spalart, 2003) have shown (and we
confirm it later) that this third mechanism is quite weak, because capturing
the overall development of the turbulent structure from a different virtual
origin is more important than applying non-zero boundary conditions for the
turbulent quantities.

###  4.2 Combining elevation and roughness variations

To examine the relative strength and the combination of these three source
mechanisms, we first perform linearised calculations for $\Rey_{\tau}=1000$
for a smooth sinusoidal wall where $f(x_{3})=\cos(2\pi/\Lambda x_{3})$, with
period $\Lambda=1$, by activating only one source mechanism each time. Results
are reported in figure 9. The wall-normal velocity profiles, taken at
$x_{3}=0$ over the crest of the topography, show that source mechanism $A$
produces a down-welling flow over the region where the slip velocity is
negative. Conversely, the decrease of the distance from the wall over the
crest, source mechanism $B$, produces an up-welling of slightly greater
magnitude, while source mechanism $B^{\prime}$ is much weaker than the first
two, as discussed.

Figure 9: Wall-normal velocity profiles over the crest of a sinusoidal
topography, for $\Lambda=1$ and $\Rey_{\tau}=1000$, obtained by activating
each of the three secondary-flow-inducing source mechanisms in turn.

Because the proposed model is linear, the superposition principle applies and
the effect of varying simultaneously the roughness properties and the
elevation can be obtained easily by combining appropriately solutions obtained
in the two cases. For ridge-type roughness, source mechanisms $A,B$ and
$B^{\prime}$ are all active, while for strip-type roughness only source
mechanisms $B$ and $B^{\prime}$ should be retained, after inverting the sign
of the induced flow given the different orientations of $f(x_{3})$ and
$d_{0}^{(1)}$ (see figure 11).

To better demonstrate the relative importance of these mechanisms, we then
consider a configuration where strips and rectangular ridges are arranged in
phase, and the high-roughness regions are placed over the ridges (see figure
10). The width $S_{l}$ coincides with the gap between the ridges while $S_{h}$
coincides the ridge width (denoted as $G$ and $W$ in Zampino et al. (2022)).
To characterise the relative strength of the two effects, we introduce the
parameter $\beta$ as the ratio between the displacement of the virtual origin
produced by the strips and the topography, so that

$d_{0}^{(1)}(x_{3})=\beta f(x_{3}),$ (25)

with the caveat that positive displacements are in different directions
depending on the type of heterogeneity. Case $\beta=0$ corresponds to smooth
ridges, leading to secondary flows produced by lateral variations of the
elevation. As discussed in Zampino et al. (2022), the linearised RANS model
predicts an upwelling over the high-elevation regions. On the other hand, case
$\beta=1$ corresponds to the combination of the two roughness heterogeneities
where the downward displacement of the virtual origin produced by the shear-
increasing roughness is, in theory, fully compensated by a increased
elevation.

Figure 10: Superposition of high and low roughness strips with smooth
rectangular ridges modelling the effect of lateral variation of elevation over
the bottom wall of the channel. The lateral variation of the virtual origin of
the rough surface is captured by $d_{0}^{(1)}(x_{3})$ (positive if the virtual
origin is below the mean height at $x_{2}=0$), while the lateral variation of
the elevation is captured by the function $f(x_{3})$ (positive if it is above
$x_{2}=0$).

Figure 11 shows results of this analysis, where we conduct calculations over
surfaces with combined roughness and elevation as a function of $\beta$. The
kinetic energy density obtained at each composite surface is normalised with
the reference value at $\beta=0$ and is shown in panel (g). Results are shown
for $\Rey_{\tau}=1000$ and $h_{s}^{(0)}=180$, for strip (and ridge) widths
configurations characterised by constant $S_{h}=S_{l}=0.7$. The kinetic energy
shows a minimum for $\beta\approx 0.3$, where the cross-stream velocity
components vanish and no secondary flows are predicted. The location of the
minimum does not seem to depend greatly on the strip configuration (not shown
for brevity). For $\beta>0.3$ the effect of the lateral variation of the
roughness becomes dominant and the associated kinetic energy density can be
several times higher than the reference value. The resulting flow structure as
$\beta$ is increased, i.e. as the effect of roughness is increased, is shown
in panels (a-d). Maps of the perturbation streamwise and wall-normal velocity
components are shown, for $\beta=0$, in the region where the ridge-type
roughness is dominant, $\beta=0.45$ close to the minimum of
$\mathcal{K}(\beta)/\mathcal{K}(0)$, and 1, where the strip-type roughness is
dominant. For $\beta=0$ (left column), the flow topology shows an upwelling
over the ridges, as the upward displacement of the wall produced the ridge-
type roughness in this region is dominant. For $\beta=0.45$ (central column),
secondary currents have changed in direction but their strength is relatively
weak. For $\beta=1$ (right column), the flow topology shows a downwelling over
the high roughness strip, as the effect of the downward displacement of the
virtual origin produced by the roughness prevails in this region over the
influence of the increased elevation.

Figure 11: Contours of the perturbation streamwise (a-c) and wall-normal (d-f)
velocity components for combinations of ridge-type and strip-type roughness
for $S_{h}=S_{l}=0.7$. The parameter $\beta$ is 0.0, 0.45 and 1, as we move
from left to right with configurations where the effect of the elevation
heterogeneity dominates to cases where the roughness heterogeneity prevails.
Normalised kinetic energy density of the combination of both ridge-type and
strip-type roughness as a function of the parameter $\beta$. The Reynolds
number is fixed at 1000 and $h_{s}^{(0)}$ = 180.

The key result of this analysis is that lateral variations of the roughness
properties and of the elevation do not have the same impact on the strength of
the resulting secondary motions, because of the “damping” effect of the
spanwise variation of the slip velocity produced by source mechanism $A$,
active for heterogeneous elevation but not for heterogeneous roughness. The
location of the minimum of $\mathcal{K}(\beta)$ suggests that the secondary-
flow-inducing effect of roughness is about three times stronger than that of
ridges, for the same displacement of the virtual origin in absolute terms.
However, this prediction is clearly no better than the prediction of the
strength of secondary flows for the two heterogeneities. The comparison with
the heterogeneous roughness simulations of Chung et al. (2018) reported
earlier suggests that the linear model captures quite well the strength of
secondary motions over such surfaces. However, for ridges, recent work (Castro
& Kim, 2024) has suggested that the linear model over-predicts the intensity
of secondary motion for tall ridges, owing to the importance of nonlinear
effects near the corners of the ridges, while predictions can be more accurate
for short ridges that do not protrude excessively in the wall-bounded flow.
Overall, this indicates that the response of the wall-bounded flow to a
perturbation of the surface elevation is far from being linear. To our
knowledge, evidence for this claim was perhaps first given in Wang & Cheng
(2006) (see their figure 20), who showed how the vertical velocity produced by
secondary motions saturates rather quickly as the height of the topography is
increased. This last piece of evidence suggests that the secondary-flow-
inducing effect of the roughness might be even stronger than what the linear
model suggest here, although a precise quantification might require dedicated
experimental work. Overall, this could also justify the recent results of
Frohnapfel et al. (2024) who considered the same surface arrangement
considered in this section, with roughness strips located in phase with the
ridges. These authors increased the ridge height but did not observe a
reversal of the flow direction above the ridges, which was dominated by the
down-welling caused by the roughness. Using their data, we calculate the
height of the ridges to be 5.16% the height of the channel and the downward
displacement of the virtual origin to be 1.08%, from their homogeneous
roughness data, for a ratio $\beta\simeq 0.21$. This ratio is clearly to the
left of the minimum in figure 11, reinforcing the idea that the response of
the wall bounded flow to finite lateral variations of the elevation are less
intense than predicted by a linear model for infinitesimal perturbations.

##  5 Discussion and concluding remarks

In this paper, we proposed a linearised-RANS-based framework to predict the
structure of Prandlt’s secondary flows of the second kind developing over
laterally heterogeneous rough surfaces. The work extends our previous efforts
on modelling smooth non-planar surfaces (Zampino et al., 2022), e.g. surfaces
with longitudinal ridges. The model couples the linearised RANS equations with
the Spalart-Allmaras (SA) transport equation to capture the altered turbulent
structure. Rough surfaces with alternating streamwise-aligned strips of high
and low roughness are modelled using established RANS modelling strategies
available in the literature for homogeneous rough surfaces (Aupoix, 2007). In
a nutshell, these strategies adopt a virtual origin framework, whereby the
shift of the logarithmic profile is obtained by displacing beneath the
boundary of the numerical domain the virtual origin (Rotta, 1962) from which
the turbulent quantities (e.g. the turbulent viscosity in the SA transport
model) develop. This results in altered boundary conditions as well as a
domain forcing term when the distance from the wall appears in the turbulence
model’s transport equations. The framework also employs a nonlinear Reynolds
stress model, i.e. the Quadratic Constitutive Relation (QCR, see Spalart
(2000)), so that secondary currents induced by the inhomogeneity of
anisotropic turbulent stresses can be predicted (Speziale, 1982).

There are several aspects that the model predicts remarkably well in agreement
with previous observations. One first aspect is the presence of three separate
flow regimes as the strip width is increased (Chung et al., 2018). For narrow
strips, the linear model support a flow structure consisting of rolls
localised in the vicinity of the wall and having wall-normal size scaling
linearly with the strip width. For wider strips, the flow structure tends to
an “isolated-vortex” regime, where streamwise vorticity is concentrated in a
roughly square region localised around the transition between strips, while
the bulk flow in the centre of the roughness strips tends to its homogeneous
rough-wall flow dictated by the local roughness properties. In the
intermediate regime, secondary currents are most intense when the high and low
roughness strips have the same widths, and are equal to about $0.7$ the half-
channel height. The model also provides adequate quantitative predictions of
the intensity of the cross-stream velocity component, compared to, e.g., the
numerical simulations of (Chung et al., 2018). A second aspect concerns
tertiary structures and flow reversal, not observed in previous studies on
roughness strips that have most often examined high and low roughness strips
of equal width. The linear model predicts that these phenomena only appear
when the strips have different width. For instance, for $\Lambda=2.8$, flow
reversal is strongest on the low-roughness strip when this strip is about 4
times wider than the high-roughness strip. It would be interesting to confirm
this prediction through experiments or simulations. A third aspect concerns
the occurrence of low- and high- momentum pathways flanking the longitudinal
rolls, where high-speed flow may be found on the high-roughness strip (and
vice-versa) in regions dominated by the vertical velocities induced by the
rolls. Away from the rolls, or for wide strips, the expected relationship
between surface roughness and streamwise velocity defect is recovered.

From a practical standpoint, the advantage of the present approach is,
undoubtedly, its computational efficiency. However, the ability to rapidly
probe the parameter space has enabled progress to be made on a more
fundamental standpoint. Specifically, previous work has suggested that
secondary currents may be the time-averaged picture of naturally-occurring
large scale motions locked in place by the heterogeneity. The robustness of
the above-mentioned similarities between the flow structure predicted by the
present framework and previous observations suggests an alternative input-
output perspective whereby these currents are the output response of the
Navier-Stokes operator linearised about the turbulent mean and subjected of a
steady, streamwise-independent forcing localised at the wall, associated to
the lateral perturbation of the surface characteristics. This perspective
complements the well accepted viewpoint that instantaneous coherent structures
in wall-bounded turbulence may be described to a satisfactory degree by the
output properties of the linearised operator subjected to a random forcing
(McKeon & Sharma, 2010; Hwang & Cossu, 2010). Admittedly, this viewpoint does
not explain the observed interplay between large scale motions and secondary
structures (Wangsawijaya & Hutchins, 2022), whereby energy of the former
leeches into the latter. One possible explanation that is worth exploring
further is that the mean flow distortion produced by the secondary currents
may locally alter the selective amplification properties of the linearised
operator, producing a spanwise modulation of the nature of large-scale motions
that may be interpreted as an energy interaction between secondary currents
and large-scale motions.

A second key output of this study is that it offers a unified perspective to
examine both ridge-type and strip-type roughness. Examination of these two
cases within the present framework shows that secondary motions produced by
complex surface heterogeneities, e.g. arbitrary combinations of elevation and
roughness properties, may be seen as originating from two separate source
mechanisms. The first is a lateral variation of the virtual origin from which
the mean turbulence structure develops. The sign of this variation is opposite
for ridge-type and strip-type roughness: the virtual origin is shifted upwards
by ridges, and downwards by higher roughness. The second source mechanism is
mediated by the lateral variation of the streamwise slip velocity in the
vicinity of the wall, and the associated spanwise gradients of the streamwise
velocity. This source mechanism captures the acceleration/deceleration
perceived by the bulk flow above the troughs and crests of a non-planar
topography, respectively, or when the mean roughness height varies laterally.
In ridge-type roughness, we have shown that this source mechanism damps the
first so that the resulting secondary motions are weaker compared to those
that would be predicted from the first mechanism. In other words, _for the
same lateral variation of the virtual origin_ strip-type roughness produces
more intense secondary flows than ridge-type roughness. The caveat is that
this perspective applies, in the limit where the lateral variation of surface
attributes is “small”, in the region of validity of the linear model. For
finite amplitude perturbations of the surface attributes, these predictions
would need to be further verified. In any case, the present modelling
framework suggests that the mean roughness height, a geometric quantity used
in previous studies that considered the coupling of roughness and elevation is
not the important quantity to be monitored when investigating combinations of
elevation and roughness. Instead, we suggest that the notion of the virtual
origin, a dynamical parameter associated to the downward shift of the
logarithmic distribution, should be considered.

Like all other analyses of the linearised Navier-Stokes operator, the present
approach yields useful insight but has its limitations. The ability of a one-
equation turbulence equipped with a nonlinear Reynolds stress model to capture
the unsteady, meandering motion of secondary structure may be questioned. A
more extensive evaluation of RANS-based methods, including turbulence models,
Reynolds stress models, compared with high-fidelity simulations is certainly
warranted, although it must be pointed out that the SA-QCR does seem to
capture fairly well the anisotropic nature of the Reynolds stress tensor in
square-duct flow (Modesti, 2020). It is also unclear how sensible is the
streamwise independence assumption used here when secondary motions display a
strong meandering (Kevin et al., 2019). Further, given that predictions are in
good agreement with experimental results, linear mechanisms may be in fact
dominant. However, the importance of nonlinear effects in altering the
structure of secondary flows is not clearly understood. Analyses of the
streamwise vorticity budget for ridge-type roughness (Castro & Kim, 2024)
shows that nonlinear effects appear significant near sharp geometric features,
and similar conclusions might apply to regions where roughness properties vary
sharply. Nonlinear effects must also be introduced in order to predict the
change in drag (Zampino, 2023), since a third limitation of the present linear
approach is that it offers no insight into the dependence of drag on the
surface properties. This would be highly useful for predicting the drag of
real-world surfaces, which is currently a topic of active research. The
perturbation velocities calculated at first order do allow dispersive stresses
and the associated increase in drag to be evaluated. However, the change in
total drag is a second order term of the difference between the roughness of
the two strips and additional terms not captured by the linear model are
required.

Data access statement: All data supporting this study will be openly available
from the University of Southampton repository.

Acknowledgments: We acknowledge funding from EPSRC (Grant ref: EP/V00199X/1).
We are grateful for the insightful discussions with Prof Ian Castro and Prof
Jae-Wook Kim.

Declaration of Interests: The authors report no conflict of interest.

## References

  * del Álamo & Jiménez (2006) del Álamo, J. C. & Jiménez, J. 2006 Linear energy amplification in turbulent channels. J. Fluid Mech. 559, 205–213.
  * Anderson et al. (2015) Anderson, W., Barros, J. M., Christensen, K. T. & Awasthi, A. 2015 Numerical and experimental study of mechanisms reponsible for turbulent secondary flows in boundary layer flows over spanwise heterogeneous roughness. J. Fluid Mech. 768, 316–347.
  * Aupoix (2007) Aupoix, B. 2007 A General Strategy to Extend Turbulence Models to Rough Surfaces: Application to Smith’s k-L Model. Journal of Fluids Engineering 129 (10), 1245–1254.
  * Aupoix & Spalart (2003) Aupoix, B. & Spalart, P. R. 2003 Extensions of the spalart-allmaras turbulence model to account for wall roughness. Int. J. Heat Fluid Flow 24, 454–462.
  * Bai et al. (2018) Bai, H. L., Kevin, Hutchins, H. & Monty, J. P. 2018 Turbulence modifications in a turbulent boundary layer over a rough wall with spanwise-alternating roughness strips. Phys. Fluid 30, 055105.
  * Barros & Christensen (2014) Barros, J. M. & Christensen, K. T. 2014 Observation of turbulent secondary flows in a rough-wall boundary layer. J. Fluid Mech. 748, R1.
  * Bottaro et al. (2006) Bottaro, A., Soueid, H. & Galletti, B. 2006 Formation of secondary vortices in turbulent square-duct flow. AIAA Journal 44, 803–811.
  * Bou-Zeid et al. (2007) Bou-Zeid, E., Parlange, M. B. & Meneveau, C. 2007 On the parameterization of surface roughness at regional scales. Journal of the atmospheric sciences 64 (1), 216–227.
  * Brandt (2014) Brandt, L. 2014 The lift-up effect: The linear mechanism behind transition and turbulence in shear flows. Eur. J. Mech. B Fluids 47, 80–96.
  * Brereton & Yuan (2018) Brereton, G.J. & Yuan, J. 2018 Wall-roughness eddy viscosity for Reynolds-averaged closures. International Journal of Heat and Fluid Flow 73, 74–81.
  * Castro & Kim (2024) Castro, I. P. & Kim, J.W 2024 Secondary motions in turbulent ribbed channel flows. J. Fluid Mech. 988, A2.
  * Castro et al. (2021) Castro, I. P, Kim, J. W., Stroh, A. & Lim, H. C. 2021 Channel flow with large longitudinal ribs. J. Fluid Mech. 915, A92.
  * Chen et al. (1997) Chen, W. L., Lien, F. S. & Leschziner, M. A. 1997 Non-linear eddy-viscosity modelling of transitional boundary layers pertinent to turbomachine aerodynamics. Int. J. Heat Fluid Flow 19, 297–306.
  * Chernyshenko & Baig (2005) Chernyshenko, S. I. & Baig, M. F. 2005 The mechanism of streak formation in near-wall turbulence. J. Fluid Mech. 544, 99–131.
  * Chung et al. (2018) Chung, D., Monty, J. P. & Hutchins, N. 2018 Similarity and structure of wall turbulence with lateral wall shear stress variations. J. Fluid Mech. 847, 591–613.
  * Cossu et al. (2009) Cossu, C., Pujals, G. & Depardon, S. 2009 Optimal transient growth and very large-scale structures in turbulent boundary layers. J. Fluid Mech. 619, 79–94.
  * Durbin et al. (2000) Durbin, P. A., Medic, G., Seo, J.-M., Eaton, J. K. & Song, S. 2000 Rough Wall Modification of Two-Layer $\kappa-\epsilon$. J. Fluids Eng. 123, 16–21.
  * Forooghi et al. (2020) Forooghi, P., Yang, X.I.A. & Abkar, M. 2020 Roughness-induced secondary flows in stably stratified turbulent boundary layers. Physics of Fluids 32 (10), 105118.
  * Frohnapfel et al. (2024) Frohnapfel, B., von Deyn, L., Yang, J., Neuhauser, J., Stroh, A., Örl”̈, R. & Gatti, D. 2024 Flow resistance over heterogeneous roughness made of spanwise-alternating sandpaper strips. J. Fluid Mech. 980, A31.
  * Goldstein & Tuan (1998) Goldstein, D. B. & Tuan, T.-C. 1998 Secondary flows induced by riblets. J. Fluid Mech. 363, 115–151.
  * Grigson (1992) Grigson, C. 1992 Drag losses of new ships caused by hull finish. Journal Ship Res. 36, 182–196.
  * Hinze (1973) Hinze, J. O. 1973 Experimental investigation on secondary currents in the turbulent flow throughout a straight conduit. App. Sci. Res. 28 (9), 453–465.
  * Hutchins et al. (2023) Hutchins, N., Ganapathisubramani, B., Schultz, M.P. & Pullin, D.I. 2023 Defining an equivalent homogeneous roughness length for turbulent boundary layers developing over patchy or heterogeneous surfaces. Ocean Engineering 271, 113454.
  * Hwang & Lee (2018) Hwang, H. G. & Lee, J. H. 2018 Secondary flows in turbulent boundary layers over longitudinal surface roughness. Phys. Rev. Fluids 3, 014608.
  * Hwang & Cossu (2010) Hwang, Y. & Cossu, C. 2010 Linear non-normal energy amplification of harmonic and stochastic forcing in the turbulent channel flow. J. Fluid Mech. 664, 51–73.
  * Jiménez (2004) Jiménez, J. 2004 Turbulent flows over rough walls. Annu. Rev. Fluid Mech. 36, 173–196.
  * Kevin et al. (2019) Kevin, Monty, J. & Hutchins, N. 2019 The meandering behaviour of large-scale structures in turbulent boundary layers. J. Fluid Mech. 865, R1.
  * Knoll & Keyes (????) Knoll, D. A. & Keyes, D. E. ???? Jacobian-free Newton-Krylov methods: a survey of approaches and applications. J. Comput. Phys 193, 357–397.
  * Lee et al. (2019) Lee, Jae Hwa, Sung, Hyung Jin & Adrian, Ronald J. 2019 Space–time formation of very-large-scale motions in turbulent pipe flow. J. Fluid Mech. 881, 1010–1047.
  * Long et al. (2023) Long, Y., Wang, J. & Pan, C. 2023 The influence of roughness-element-spacing on turbulent entrainment over spanwise heterogeneous roughness. Physics of Fluids 35 (8), 085127\.
  * McKeon & Sharma (2010) McKeon, B. J. & Sharma, A. S. 2010 A critical-layer framework for turbulent pipe flow. J. Fluid Mech. 658, 336–382.
  * Medjnoun et al. (2018) Medjnoun, T., Vanderwel, C. & Ganapathisubramani, B. 2018 Characteristics of turbulent boundary layers over smooth surfaces with spanwise heterogeneities. J. Fluid Mech. 838.
  * Medjnoun et al. (2020) Medjnoun, T., Vanderwel, C. & Ganapathisubramani, B. 2020 Effects of heterogeneous surface geometry on secondary flows in turbulent boundary layers. J. Fluid Mech. 886, A31.
  * Mejia-Alvarez & Christensen (2013) Mejia-Alvarez, R. & Christensen, K. 2013 Wall-parallel stereo particle-image velocimetry measurements in the roughness sublayer of turbulent flow overlying highly irregular roughness. Phy. Fluids 25, 015106.
  * Meyers et al. (2019) Meyers, J., Ganapathisubramani, B. & Cal, R. B. 2019 On the decay of dispersive motions in the outer region of rough-wall boundary layers. J. Fluid Mech. 862, R5.
  * Modesti (2020) Modesti, D 2020 A priori tests of eddy viscosity models in square duct flow. Theoretical and Computational Fluid Dynamics 34 (5), 713–734.
  * Morra et al. (2019) Morra, P., Semeraro, O., Henningson, D. S. & Cossu, C. 2019 On the relevance of reynolds stresses in resolvent analyses of turbulent wall-bounded flows. Journal of Fluid Mechanics 867, 969–984.
  * Neuhauser et al. (2022) Neuhauser, J., Schäfer, K., Gatti, D. & Frohnapfel, B. 2022 Simulation of turbulent flow over roughness strips. J. Fluid Mech. 945, A14.
  * Perkins (1970) Perkins, H. J. 1970 The formation of streamwise vorticity in turbulent flow. J. Fluid Mech. 44, 721–740.
  * Pickering et al. (2021) Pickering, Ethan, Rigas, Georgios, Schmidt, Oliver T., Sipp, Denis & Colonius, Tim 2021 Optimal eddy viscosity for resolvent-based models of coherent structures in turbulent jets. Journal of Fluid Mechanics 917, A29.
  * Prakash & Laurendeau (2020) Prakash, A. & Laurendeau, E. 2020 Consistent surface roughness extension for wall functions. Int. J. Heat Fluid Flow 82, 108552.
  * Prandtl (1952) Prandtl, L. 1952 Essentials of fluid dynamics. Hafner.
  * Pujals et al. (2009) Pujals, G., Garcìa-Villalba, M., Cossu, C. & Depardon, S. 2009 A note on optimal transient growth in turbulent channel flows. Phy. Fluids 21, 015109.
  * Reynolds & Hussain (1972) Reynolds, W. C. & Hussain, A. K. M. F. 1972 The mechanism of an organized wave in turbulent shear flow. part 3. theoretical models and comparisons with experiments. J. Fluid Mech. 54, 263–288.
  * Rotta (1962) Rotta, J. C. 1962 Turbulent boundary layers in incompressible flow. Progress in Aerospace Sciences 2, 1–95.
  * Schäfer et al. (2022) Schäfer, K., Stroh, A., Forooghi, P. & Frohnapfel, B. 2022 Modelling spanwise heterogeneous roughness through a parametric forcing approach. J. Fluid Mech. 930, A7.
  * Spalart (2000) Spalart, P. R. 2000 Strategies for turbulence modelling and simulations. Int. J. Heat Fluid Flow 21, 252–263.
  * Spalart & Allmaras (1994) Spalart, P. R. & Allmaras, S. R. 1994 A one-equation turbulence model for aerodynamic flows. Rech. Aerosp. 1, 5–21.
  * Speziale (1982) Speziale, C. G. 1982 On turbulent secondary flows in pipes of noncircular cross-section. Int. J. Engng Sci. 20 (7), 863–872.
  * Speziale et al. (1991) Speziale, C. G., Sarkar, S. & Gatski, T. B. 1991 Modelling the pressure-strain correlation of turbulence: an invariant dynamical system approach. J. Fluid Mech. 227, 254–272.
  * Stroh et al. (2020) Stroh, A., Schäfer, K., Frohnapfel, B. & Forooghi 2020 Rearrangement of secondary flow over spanwise heterogenous roughness. J. Fluid Mech. 885, R5.
  * Stroh et al. (2016) Stroh, A., Y. Hasegawa, J. Kriegseis & Frohnapfel, B. 2016 Secondary vortices over surfaces with spanwise varying drag. J. Turbul. 17, 1142–1158.
  * Suga et al. (2006) Suga, K., Craft, T. J. & Iacovides, H. 2006 An analytical wall-function for turbulent flows and heat transfer over rough walls. Int J Heat Fluid Flow 27 (5), 852–866.
  * Townsend (1976) Townsend, A. A. R. 1976 The structure of turbulent shear flow. Cambridge university press.
  * Vanderwel & Ganapathisubramani (2015) Vanderwel, C. & Ganapathisubramani, B. 2015 Effects of spanwise spacing on large-scale secondary flows in rough-wall turbulent boundary layers. J. Fluid Mech. 774, R2.
  * Vidal et al. (2018) Vidal, A., Nagib, H. M., Schlatter, P. & Vinuesa, R. 2018 Secondary flow in spanwise-periosic in-phase sinusoidal channels. J. Fluid Mech. 851, 288–316.
  * Viswanath (2007) Viswanath, D. 2007 Recurrent motions within plane couette turbulence. J. Fluid Mech. 580, 339–358.
  * Wang & Cheng (2006) Wang, Z. Q. & Cheng, N. S. 2006 Time-mean structure of secondary flows in open channel with longitudinal bedforms. Advances in Water Resources 29, 1634–1649.
  * Wangsawijaya & Hutchins (2022) Wangsawijaya, D.D. & Hutchins, N. 2022 Investigation of unsteady secondary flows and large-scale turbulence in heterogeneous turbulent boundary layers. J. Fluid Mech. 934, A40.
  * Wangsawijaya et al. (2020) Wangsawijaya, D. D., Baidya, R., Chung, D., Marusic, I. & Hutchins, N. 2020 The effect of spanwise wavelength of surface heterogeneity on turbulent secondary flows. J. Fluid Mech. 894, A7.
  * Willingham et al. (2014) Willingham, D., Anderson, W., Christensen, K. T. & Barros, J. M. 2014 Turbulent boundary layer flow over transverse aerodynamic roughness transitions: induced mixing and flow characterization. Phy. Fluids 26, 025111.
  * Yang & Anderson (2018) Yang, J. & Anderson, W. 2018 Numerical study of turbulent channel flow over surfaces with variable spanwise heterogeneities: Topographically-driven secondary flows affect outer-layer similarity of turbulent length scales. Flow Turb. Comb. 100, 1–17.
  * Zampino (2023) Zampino, G. 2023 Linearized reynolds-averaged navier-stokes model for the prediction of secondary flows on heterogeneous surfaces. PhD thesis, University of Southampton.
  * Zampino et al. (2022) Zampino, G., Lasagna, D. & Ganapathisubramani, B. 2022 Linearised Reynolds-averaged predictions of secondary currents in turbulent channels with topographic heterogeneity. J. Fluid Mec. 944, A4.
  * Zampino et al. (2023) Zampino, G., Lasagna, D. & Ganapathisubramani, B. 2023 Scaling of secondary flows with surface parameters: A linear approach. Int. J. Heat and Fluid Flow 102, 109148.
  * Zampiron et al. (2020) Zampiron, A., Cameron, S. & Nikora, V. 2020 Secondary currents and very-large-scale motions in open-channel flow over streamwise ridges. J. Fluid Mech. 887, A17.
  * Zampiron et al. (2021) Zampiron, A., Cameron, S. & Nikora, V. 2021 Momentum and energy transfer in open-channel flow over streamwise ridges. J. Fluid Mech. 915, A42.
  * Zhdanov et al. (2024) Zhdanov, O., Jelly, T. O & Busse, A. 2024 Influence of ridge spacing, ridge width, and reynolds number on secondary currents in turbulent channel flow over triangular ridges. Flow, Turbulence and Combustion 112 (1), 105–128.