1. Introduction
Scalar transport and transfer near the interfaces between a turbulent flow and a porous medium are crucial in a variety of natural systems and engineering applications. Among the technical applications are fuel cells (e.g. Najmi et al. Reference Najmi, El-Tabach, Chetehouna, Gascoin and Falempin2016; Yang et al. Reference Yang, Xu, Li, Chen, Hu, He and Du2023), cooling processes (e.g. Dahmen et al. Reference Dahmen, Müller, Rom, Schweikert, Selzer and von Wolfersdorf2015; Ranjbar, Pouransari & Siavashi Reference Ranjbar, Pouransari and Siavashi2021) and drying processes (e.g. Sharma et al. Reference Sharma, Chen and Vu Lan2009; Zhao et al. Reference Zhao, Rui, Rong and Wang2021). In natural systems, notable examples include the evaporation from soils (e.g. Mosthaf et al. Reference Mosthaf, Baber, Flemisch, Helmig, Leijnse, Rybak and Wohlmuth2011; Kiemle et al. Reference Kiemle, Heck, Coltman and Helmig2023) as well as processes within the hyporheic zone of a porous river bed (e.g. Lewandowski et al. Reference Lewandowski2019; Tewari, Singh & Gaur Reference Tewari, Singh and Gaur2021). The hyporheic zone plays a decisive role in the context of river ecosystems. It encompasses parts of the sediment bed where the interaction between turbulent streamflow and pore water facilitates a bidirectional exchange of energy, momentum and mass (e.g. Boano et al. Reference Boano, Harvey, Marion, Packman, Revelli, Ridolfi and Wörman2014; Woessner Reference Woessner, Hauer and Lamberti2017). Hyporheic transport processes at the grain scale are primarily driven by hydrodynamic forces (Boano et al. Reference Boano, Harvey, Marion, Packman, Revelli, Ridolfi and Wörman2014). These processes are significant due to their occurrence over large areas and their impact within regions of high biogeochemical activity directly below the sediment–water interface (Marion et al. Reference Marion2014). The sediment–water interface can be structured in several layers (e.g. Nikora et al. Reference Nikora, Goring, McEwan and Griffiths2001), in which different scalar transport processes are significant. In the water column far above the sediment bed, turbulent mixing is dominant, while at the interface and below, dispersive and diffusive scalar transport gain importance (Voermans, Ghisalberti & Ivey Reference Voermans, Ghisalberti and Ivey2018). In each layer, different scalings are appropriate, which makes a general description of transport across the interface a challenging task.
The Reynolds analogy (Reynolds Reference Reynolds1874, Reference Reynolds1961) can serve as a basis to relate scalar transport to momentum transport, as it implies that both quantities are transferred by the same mechanisms. Kim & Moin (Reference Kim and Moin1989) were among the first to investigate passive scalar transport at
${\textit{Pr}} \leqslant 2$
using direct numerical simulation (DNS). The scalar was introduced via a volume source and absorbed on the smooth and impermeable channel walls, which facilitated a direct comparison between velocity und scalar field. In the near-wall region, a strong correlation between streamwise velocity fluctuations and scalar fluctuations was observed, with both exhibiting a similar pattern of streamwise elongated streaks. Subsequent studies of scalar transport in smooth-wall channels targeted higher Reynolds numbers (e.g. Abe, Kawamura & Matsuo Reference Abe, Kawamura and Matsuo2004; Pirozzoli, Bernardini & Orlandi Reference Pirozzoli, Bernardini and Orlandi2016), whereas others expanded the parameter space towards higher Schmidt or Prandtl numbers for the scalar transport (e.g. Kawamura et al. Reference Kawamura, Ohsaka, Abe and Yamamoto1998; Schwertfirm & Manhart Reference Schwertfirm and Manhart2007). For the range of
${\textit{Pr}} \in [1,10]$
, Alcántara-Á vila & Hoyas (Reference Alcántara-Á vila and Hoyas2021) confirmed that the turbulent diffusivity and turbulent viscosity have similar values in the far-wall region.
As summarised in the reviews of Hantsis & Piomelli (Reference Hantsis and Piomelli2024) and Kadivar & Garg (Reference Kadivar and Garg2025), numerous studies have addressed passive scalar transport over rough, yet impermeable, surfaces. Most numerical studies consider Schmidt numbers
${\textit{Sc}} = \nu \kern-0.5pt/\kern-1pt \varGamma$
or Prandtl numbers
${\textit{Pr}} = \nu \kern-0.5pt/\kern-1pt \alpha$
near unity, which implies that the smallest scales in the scalar field (Batchelor scale) are of the same size as the smallest scales of turbulent motion (Kolmogorov scale). Here,
$\nu$
is the kinematic viscosity, while
$\varGamma$
and
$\alpha$
represent the molecular or thermal diffusivity, respectively. From the underlying equations, it can be expected that the Reynolds analogy breaks down if pressure gradients play a decisive role. Several studies confirmed that the Reynolds analogy is not valid in the proximity of sufficiently rough surfaces (e.g. Dipprey & Sabersky Reference Dipprey and Sabersky1963; Leonardi et al. Reference Leonardi, Orlandi, Djenidi and Antonia2015; Li & Bou-Zeid Reference Li and Bou-Zeid2019; MacDonald, Hutchins & Chung Reference MacDonald, Hutchins and Chung2019; Kuwata Reference Kuwata2021). To quantify the progressive breakdown of the analogy with increasing roughness Reynolds number
$k^+_s$
(the inner-scaled equivalent sand-grain roughness height), Forooghi, Stripf & Frohnapfel (Reference Forooghi, Stripf and Frohnapfel2018) resorted to the Reynolds analogy factor. In addition to the Reynolds analogy factor, the turbulent Schmidt number can be employed to assess the analogy between momentum and scalar transfer. Slightly different definitions of this number appear in the literature (e.g. Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2016; MacDonald et al. Reference MacDonald, Hutchins and Chung2019; Hantsis & Piomelli Reference Hantsis and Piomelli2024). By including the dispersive transport in addition to the turbulent transport, an effective Schmidt number can be obtained (e.g. Kuwata Reference Kuwata2021; Hantsis & Piomelli Reference Hantsis and Piomelli2022). Following the argumentation of Abe, Antonia & Kawamura (Reference Abe, Antonia and Kawamura2009), Secchi et al. (Reference Secchi, Gatti, Piomelli and Frohnapfel2025) recently suggested a novel framework for the investigation of the Reynolds analogy based on budget terms. In the turbulent free-flow region above the roughness sublayer, rough-wall scalar statistics were found to show similarities with their smooth-wall equivalents (e.g. MacDonald et al. Reference MacDonald, Hutchins and Chung2019; Peeters & Sandham Reference Peeters and Sandham2019; Hantsis & Piomelli Reference Hantsis and Piomelli2020, Reference Hantsis and Piomelli2022). Based on these observations, Hantsis & Piomelli (Reference Hantsis and Piomelli2024) concluded in their review that Townsend’s similarity hypothesis is valid also for scalar transport at sufficiently high Reynolds numbers, whereas Kadivar & Garg (Reference Kadivar and Garg2025) pointed out discrepancies between smooth- and rough-wall scalar profiles.
In contrast to the large body of rough-wall studies, high-resolution investigations of scalar transfer across porous, permeable walls are relatively rare. Chandesris et al. (Reference Chandesris, D’Hueppe, Mathieu, Jamet and Goyeau2013) investigated the effect of boundary conditions on the transport of a passive scalar in a three-dimensional Cartesian grid of cubes underneath a turbulent channel flow. They replicated the flow case of Breugem & Boersma (Reference Breugem and Boersma2005) and used a Prandtl number of
${\textit{Pr}} = 0.1$
. Comparing cases with zero and constant scalar flux through the cube surfaces (i.e. Neumann-type boundary conditions), they found that the effective turbulent diffusivity of the scalar within the porous medium remained unchanged. Nishiyama, Kuwata & Suga (Reference Nishiyama, Kuwata and Suga2020) used DNS to investigate the transport of a passive temperature field at
${\textit{Pr}}=0.71$
within a turbulent channel flow over regularly structured anisotropic porous media with high porosities ranging from
$0.6$
to
$0.8$
. Different scalar boundary conditions were applied on the surfaces of the porous media, one of them being a Dirichlet-type isothermal condition, which prescribed a fixed uniform temperature rather than a fixed scalar flux. Nishiyama et al. (Reference Nishiyama, Kuwata and Suga2020) observed that the fluid temperature adapts to the fixed surface temperature shortly below the surface of the porous medium. Compared with conjugate heat transfer between fluid and solid, the fixed surface temperature was found to reduce the turbulent and dispersive scalar transport below the surface of the porous media. This suggests that the type of scalar boundary condition has a significant influence on the scalar transport behaviour.
To investigate hyporheic exchange processes, several experimental studies (e.g. Richardson & Parr Reference Richardson and Parr1988; Nagaoka & Ohgaki Reference Nagaoka and Ohgaki1990; Packman, Salehin & Zaramella Reference Packman, Salehin and Zaramella2004; Chandler et al. Reference Chandler, Guymer, Pearson and van Egmond2016) have employed conservative tracers, which are neither involved in chemical reactions nor in sorption processes on the fluid–solid interfaces. The experimental data entered several meta-studies (e.g. O’Connor & Hondzo Reference O’Connor and Hondzo2008; Grant, Stewardson & Marusic Reference Grant, Stewardson and Marusic2012; Voermans et al. Reference Voermans, Ghisalberti and Ivey2018), which proposed effective diffusivity models. v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2025a
) simulated hyporheic scalar transport by a pore-resolved numerical simulation of turbulent flow over random sphere packs, which abstractly represent the sediment bed. In contrast to the porous media used by Chandesris et al. (Reference Chandesris, D’Hueppe, Mathieu, Jamet and Goyeau2013) or Nishiyama et al. (Reference Nishiyama, Kuwata and Suga2020), the monodisperse random sphere packs were denser and had an irregular surface. At
${\textit{Sc}} = 1$
, a bed-normal scalar flux was induced by different scalar concentrations prescribed at the bottom and at the top domain boundary. The absence of volume scalar sources combined with zero-flux scalar boundary conditions on the sphere surfaces caused the scalar to behave like a tracer. A key observation in this previous study was the strong spatial heterogeneity of the scalar field within the random sphere packs, which provided favourable conditions for dispersive scalar transport below the sediment–water interface and emphasised the influence of even small scalar-induced buoyancy effects (v. Wenczowski et al. Reference v. Wenczowski, Sakai, Casas-Mulet, Chiogna, Geist and Manhart2025). The spatial heterogeneity of the scalar is even more remarkable as
${\textit{Sc}} = 1$
implies a higher scalar diffusivity than a Schmid number of
${\textit{Sc}} = O(10^3)$
, which is commonly assumed for the transport of dissolved substances in water.
The above-mentioned observations raise the question how the scalar boundary conditions on the surfaces influence the transport through the sediment and whether results obtained with conservative tracers and adiabatic boundary conditions are transferable to processes like the sorption of solutes (e.g. Gandy, Smith & Jarvis Reference Gandy, Smith and Jarvis2007), the heat transfer between solid and fluid phase (e.g. Wu et al. Reference Wu, Gomez-Velez, Krause, Singh, Wörman and Lewandowski2020) or the intake and release of substances by biofilms (e.g. Battin et al. Reference Battin, Besemer, Bengtsson, Romani and Packmann2016). The above-mentioned processes are characterised by multiple and possibly unknown parameters, such that we deliberately choose another idealised boundary condition for the present study. While zero-flux boundary conditions prevent the sorption of solutes, a Dirichlet boundary condition allows to set a uniform fixed scalar value on the sphere surfaces, while it does not restrict the absorption of scalars. Thus, all the above-mentioned processes would lie somewhere in between those two diametrically different boundary conditions. Henceforth, we will use the terms ‘isothermal’ or ‘isoscalar’ for the uniform fixed-value scalar boundary. Similarly, the term ‘adiabatic’ is used for zero-flux Neumann-type boundary conditions.
In the scope of this study, we address the following research questions: (i) Do momentum and scalar absorption exhibit similar behaviour in the near-interface region? (ii) What are the main differences compared with zero-flux scalar boundary conditions on the porous medium? (iii) Can parallels be drawn between the velocity and scalar fields in the near-interface region? Due to its focus on the previous points, the present study complements the findings reported by v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024, Reference v. Wenczowski and Manhart2025a ). The underlying equations and their formulation within the time–space double-averaging framework are found in § 2. The case configuration, the parameter space and the applied numerical methods are introduced in § 3. In § 4, we present the main findings and include comparisons to previous studies where appropriate. A discussion follows in § 5 before § 6 concludes the study.
2. Theory
In the following, the underlying equations are summarised (§ 2.1) and the double-averaging analysis framework is outlined (§ 2.2). A formulation of the double-averaged scalar transport equation is provided (§ 2.3), which distinguishes scalar fluxes within the fluid and fluxes across the fluid–solid interface. Also, we briefly discuss the budget equation for temporal scalar fluctuations (§ 2.4).
2.1. Governing equations
In the present study, the incompressible Navier–Stokes equations are solved for a Newtonian fluid with a kinematic viscosity
$\nu$
and a density
$\rho$
. The conservation of momentum and mass are formulated in (2.1) using the Einstein summation notation,
\begin{equation} \frac {\partial u_i}{\partial t} + u_{\!j} \, \frac {\partial u_i}{\partial x_{\!j}} = - \frac {1}{\rho } \frac {\partial p}{\partial x_i} + \nu \frac { \partial ^2 u_i }{ \partial x_{\!j} \partial x_{\!j} } + g_i , \quad \frac {\partial u_i}{\partial x_i} = 0. \end{equation}
The coordinates
$x_1 \equiv x$
and
$x_2 \equiv y$
refer to the streamwise and spanwise direction, respectively. The coordinate
$x_3 \equiv z$
specifies the vertical position above the sediment bed. The components of the velocity vector are denoted as
$u_1 \equiv u$
,
$u_2 \equiv v$
and
$u_3 \equiv w$
, while
$p$
is the pressure. Further, the volume force
$g_i$
acts on the fluid. The computed divergence-free flow field enters the solution of the advection–diffusion equation for a passive scalar, i.e.
\begin{equation} \frac {\partial c}{\partial t} + u_{\!j} \: \frac {\partial c}{\partial x_{\!j}} = \varGamma _c \: \frac {\partial ^2 c}{\partial x_{\!j} x_{\!j}} . \end{equation}
In (2.2),
$c$
represents the concentration of the passive scalar, while
$\varGamma _c$
is the molecular scalar diffusion coefficient. Volume sources of the scalar shall not be considered. In its dimensionless form, the advection–diffusion equation can be denoted as
\begin{equation} \frac {\partial \hat {c}}{\partial \hat {t}} + \hat {u}_{j} \: \frac {\partial \hat {c}}{\partial \hat {x}_{j}} = \frac {1}{\textit{Re} \, \textit{Sc}} \, \frac {\partial ^2 \hat {c}}{\partial \hat {x}_{j}^2} . \end{equation}
In addition to the Reynolds number
${\textit{Re}}$
, (2.3) includes the dimensionless Schmidt number
${\textit{Sc}} = \nu \kern-0.5pt/\kern-1pt \varGamma _c$
. The product of both yields the Péclet number
${\textit{Pe}} = \textit{Re} \, \textit{Sc}$
. If the passive scalar
$c$
is interpreted as a temperature field without buoyancy effects, the Prandtl number
${\textit{Pr}}$
replaces the Schmidt number
${\textit{Sc}}$
. Accordingly,
$\varGamma _c$
would then represent the thermal diffusivity.
2.2. Analysis framework
Hyporheic transport processes are commonly analysed within a temporal and spatial double-averaging framework (e.g. Giménez-Curto & Lera Reference Giménez-Curto and Lera1996; Nikora et al. Reference Nikora, Goring, McEwan and Griffiths2001; Mignot, Barthelemy & Hurther Reference Mignot, Barthelemy and Hurther2009). An arbitrary quantity
$\varphi$
undergoes a Reynolds decomposition into a temporal mean value
$\overline {\varphi }$
and a temporal fluctuation
$\varphi '$
, which is expressed as
\begin{equation} \varphi {( \boldsymbol{x},t)} = \overline {\varphi }{(\boldsymbol{x})} + \varphi ^\prime {(\boldsymbol{x},t)}, \hspace {1cm} \text{where} \hspace {1cm} {\varphi }{(\boldsymbol{x})} = \frac {1}{T} \int _{0}^{T} \varphi {(\boldsymbol{x},t)} \, {\rm d}t . \end{equation}
The temporal mean
$\overline {\varphi }$
from (2.4) is further decomposed with respect to its distribution in space. As a result,
$\langle \overline {\varphi } \rangle$
represents the intrinsic average within a horizontal plane and
$\widetilde {\overline {\varphi }}$
represents a spatial variation, defined as a deviation from the in-plane average.
\begin{equation} \overline {\varphi }{(\boldsymbol{x})} = \langle \overline {\varphi } \rangle {(z)} + \widetilde {\overline {\varphi }}{(\boldsymbol{x})} , \quad \text{where} \qquad \langle \overline {\varphi } \rangle {(z)} = \frac {1}{A_{\!f}} \iint _{A_{\!f}} \overline {\varphi }{(\boldsymbol{x})} \, {\rm d}x \, {\rm d}y .\end{equation}
According to (2.5), the intrinsic average
$\langle \overline { \varphi } \rangle$
is defined as a spatial average over the fluid-filled area
$A_{\!f}$
of the horizontal averaging plane at height
$z$
, whereas the superficial average
$\langle \overline { \varphi } \rangle ^s$
represents a spatial average over the complete base area
$A_0$
of the averaging plane. As expressed by (2.6), the in-plane porosity
$\theta {(z)}$
connects both spatial averages via
\begin{equation} \left \langle \overline {\varphi } \right \rangle ^s{(z)} = \theta {(z)} \: \left \langle \overline {\varphi } \right \rangle {(z)} \quad \text{with} \qquad \theta {(z)} = {A_{\!f}(z)} \kern-0.5pt/\kern-1pt {A_0} . \end{equation}
For an averaging plane at any vertical position
$z$
below the crests of the sphere pack, the fluid-filled area
$A_{\!f}$
is interrupted by solid-occupied areas. As a consequence, spatial derivatives and horizontal averaging do not commute and special rules apply for spatial gradients (e.g. Giménez-Curto & Lera Reference Giménez-Curto and Lera1996),
\begin{align} \displaystyle \left \langle \boldsymbol{\nabla }\overline {\varphi } \right \rangle = \left ( \begin{matrix} \displaystyle \left \langle \frac {\partial \overline {\varphi }}{\partial x} \right \rangle \\[3ex] \displaystyle \left \langle \frac {\partial \overline {\varphi }}{\partial y} \right \rangle \\[3ex] \displaystyle \left \langle \frac {\partial \overline {\varphi }}{\partial z} \right \rangle \end{matrix} \right ) = \left ( \begin{matrix} \displaystyle \frac {1}{A_{\!f}} \oint _s \: \overline {\varphi } \: \frac {n_x}{\sqrt {n_x^2 + n_y^2}} \, {\rm d}s \\[2ex] \displaystyle \frac {1}{A_{\!f}} \oint _s \: \overline {\varphi } \: \frac {n_y}{\sqrt {n_x^2 + n_y^2}} \, {\rm d}s \\[2ex] \displaystyle \frac {1}{\theta } \frac {\partial \theta \!\left \langle \overline {\varphi } \right \rangle }{\partial z} + \frac {1}{A_{\!f}} \oint _s \:\overline {\varphi } \: \frac {n_z}{\sqrt {n_x^2 + n_y^2}} \, {\rm d}s \end{matrix} \right ) \triangleq \left ( \begin{matrix} \displaystyle {\textit{BT}}_1( \overline {\varphi } ) \\[3ex] \displaystyle {\textit{BT}}_2( \overline {\varphi } ) \\[3ex] \displaystyle \frac {1}{\theta } \frac {\partial \theta \!\left \langle \overline {\varphi } \right \rangle }{\partial z} +{\textit{BT}}_3( \overline {\varphi } ) \end{matrix} \right )\! .\end{align}
The curve
$s$
in (2.7) represents the intersection of the horizontal averaging plane with solid surfaces of the sphere pack. The vector
$\boldsymbol{n} = (n_x, n_y, n_z)^T$
is the unit normal on the solid sphere surfaces and points out of the fluid-filled volume. As previously done by v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024), we refer to the curve integrals as boundary terms
${\textit{BT}}_i(\overline {\varphi })$
, while
${\textit{BT}}^{\,s}_i(\overline {\varphi }) = \theta {\textit{BT}}_i(\overline {\varphi })$
. The boundary terms account for the physical boundary conditions on fluid–solid interfaces, which are generally not met by the individual quantities
$\langle \overline {\varphi } \rangle$
and
$\widetilde {\overline {\varphi }}$
.
2.3. Double-averaged advection–diffusion equation
The advection–diffusion equation for a passive scalar (2.2) can be formulated in its conservative form, which emphasises that both advective and diffusive transport contribute to a directed scalar flux
${J}_{\!j} = u_{\!j} c - \varGamma _c \, {\partial {c}}/{\partial x_{\!j}}$
, while the divergence of the flux determines the local rate of change via
$\partial c \kern-0.5pt/\kern-1pt \partial t = -\partial J_{j} \kern-0.5pt/\kern-1pt \partial x_{\!j}$
. Applying the temporal and superficial spatial averaging operators together with the rules for spatial derivatives (2.7) yields the double-averaged form of this equation,
\begin{equation} \theta \frac {\partial \langle \overline {c} \rangle }{\partial t} \: = \: - \theta \!\left \langle \frac {\partial }{\partial x_{\!j}} \, \overline {J}_{\!j} \right \rangle \: = \: - \frac {\partial }{\partial z} \left ( \theta \!\left \langle \overline {J}_z \right \rangle \right ) \, - \, {\textit{BT}}_{\!j}^{\,s}\left ( \overline {J}_i \right )\! . \end{equation}
The formulation in (2.8) allows a distinction between the double-averaged vertical scalar flux within the fluid-filled volume
$\theta \langle \overline {J}_z \rangle$
and the flux across the fluid–solid interface, represented by the boundary term
${\textit{BT}}^{\,s}_{\!j}(\overline {J}_{\!j})$
. In the absence of a net volume flux in bed-normal direction, i.e.
$\langle \overline {w} \rangle = 0$
, the vertical scalar flux within the fluid-filled volume can be decomposed as
\begin{equation} \theta \!\left \langle \overline {J}_z \right \rangle \, = \; \underbrace { \theta \!\left \langle \overline {w^\prime c^\prime } \right \rangle }_{\text{turb.}} \; + \underbrace { \; \theta \!\left \langle \widetilde {\overline {w}} \: \widetilde {\overline {c}} \right \rangle }_{\text{disp.}} \; - \underbrace { \; \theta \!\left \langle \varGamma _c \frac {\partial \overline {c}}{\partial z} \right \rangle }_{\text{diff.}}\!. \end{equation}
Equation (2.9) allows a distinction among turbulent scalar transport, dispersive scalar transport and scalar transport due to molecular diffusion. Using the definition of the boundary term, the flux across the fluid–solid interface can be expressed by
\begin{align} {BT}^s_{\!j}( \overline {J}_{\!j} ) \; & = \; {BT}^s_{\!j}\left ( \overline {u_{\!j} c} \right ) \; + \; {BT}^s_{\!j}\left ( -\varGamma _c \: {\partial \overline {c}} \kern-0.5pt/\kern-1pt {\partial x_{\!j}} \right ) \nonumber\\[1ex] & = \; 0 \; + \; \frac {1}{A_0} \oint _{s} \, - \, \varGamma _c \frac {\partial \overline {c}}{\partial x_{\!j}} \, \frac {n_{\!j}}{\sqrt {n_x^2 + n_y^2}} \, {\rm d}s . \end{align}
Assuming that no-slip boundary conditions enforce a velocity
$u_{\!j}=0$
on the fluid–solid interface, the advective boundary term
${\textit{BT}}^{\,s}_{\!j} ( \overline {u_{\!j} c} )$
in (2.10) has zero-value. As a result, only scalar fluxes due to molecular diffusion through the solid–fluid interface are possible, which agrees with the formulation of the wall heat transfer term by Kuwata, Tsuda & Suga (Reference Kuwata, Tsuda and Suga2020). Note that other formulations are also found in the literature, which yield a different decomposition of the diffusive scalar transport (e.g. Hantsis & Piomelli Reference Hantsis and Piomelli2024). We favour the decomposition given by (2.9) and (2.10), as it allows a clear interpretation on the level of scalar fluxes.
2.4. Budget of scalar fluctuations
The double-averaged budget equation for temporal fluctuations
$\langle \overline {c' c'} \rangle \kern-1pt/2$
in the scalar field contains terms associated with production, dissipation and transport,
\begin{align} \frac {\partial }{\partial t} \left \langle \frac {1}{2} \overline {c' c'} \right \rangle &= \underbrace { { - \left \langle \overline {w' c'} \right \rangle \frac { \partial \langle \overline {c} \rangle }{\partial z} } \: { - \left \langle \overline {u_{\!j}^\prime c^\prime } \frac { \partial \widetilde {\overline {c}} }{ \partial x_{\!j} } \right \rangle } }_{\text{production ($\varPi _{\textit{grad}} + \varPi _{\textit{form}}$)} } \: \underbrace { - \: \varGamma _c \left \langle \overline { \frac {\partial c'}{\partial x_{\!j}} \frac {\partial c'}{\partial x_{\!j}} } \right \rangle }_{\text{dissipation ($\varepsilon _c$)}} \nonumber\\ &\quad - \frac {1}{\theta } \frac {\partial }{\partial z} \left ( \underbrace { \frac {\theta \!\left \langle \widetilde {\overline {w}} \: \overline {c' c'} \right \rangle }{2} }_{\text{advective ($T_{\textit{adv}}$)}} \: \underbrace { + \: \frac {\theta \!\left \langle \overline {w' c' c' } \right \rangle }{2} }_{\text{turbulent ($T_{\textit{turb}}$)}} \: \underbrace { + \: \frac {\theta \!\left \langle -\varGamma _c \frac { \partial \overline {c' c' } }{ \partial z } \right \rangle }{2} }_{\text{diffusive ($T_{\textit{diff}}$)}} \right )\! .\end{align}
As annotated in (2.11), two terms are involved in the production of temporal scalar fluctuations. The term
$\varPi _{\textit{grad}}$
is referred to as gradient production, as it represents the action of the bed-normal turbulent scalar transport against the vertical gradient of the double-averaged scalar field. The second production term
$\varPi _{\textit{form}}$
represents the form-induced production, as it can only occur if geometrical constraints introduce spatial variations in the scalar concentration field. Form-induced production results from the action of turbulent scalar transport against the gradients due to in-plane scalar variations, which are homogenised in the process. The term
$\varepsilon _c$
represents the removal of temporal scalar fluctuations due to diffusion and is, thus, equivalent to a diffusive dissipation term. The remaining terms represent advective, turbulent and diffusive transport.
3. Methods
The configuration of the simulation case (§ 3.1) and the properties of the sphere pack (§ 3.2) are outlined. The parameter space is presented (§ 3.3). A drag-based interface definition is discussed (§ 3.4) and the employed numerical methods are introduced (§ 3.5). For details on some of the above-mentioned aspects, we refer to v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024), where a convergence study and a validation for the flow field are documented, while a convergence study for the scalar field is reported by v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2025a ).
3.1. Case configuration
Turbulent open-channel flow over a mono-disperse random sphere pack is studied using grain-resolved single-domain DNS (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024). All temporal and spatial scales are fully resolved including the flow in the pore space. The overall set-up is analogous to free-surface river flow over an immobile gravel bed, where the individual spheres represent sediment grains. As illustrated in figure 1, the top boundary of the domain is a rigid lid with a free-slip boundary condition, which approximates the free surface of the water. A constant volume force
$g_x \gt 0$
drives the flow. The domain is periodic in both the streamwise
$x$
-direction and the lateral
$y$
-direction. In the statistically steady state, the boundary layer is fully developed and its thickness
$\delta$
is equal to the flow depth
$h$
. The free-slip bottom boundary of the domain intersects the sphere pack ensuring that streamwise momentum is absorbed solely by the no-slip sphere surfaces and not by the domain boundary.
Non-proportional sketch of the case configuration. A fixed scalar concentration
$c_{\textit{top}}$
(red) is specified at the free-slip lid. Within the turbulent free-flow region, the downward-oriented scalar flux (pink dashed lines) is constant. Scalar absorption happens on the sphere surfaces of the close-packed sphere pack, where a fixed concentration
$c_{\textit{pm}}$
(blue) is prescribed. All spheres have a diameter
$D$
and
$h$
represents the flow depth.
A Dirichlet-type boundary condition is applied at the rigid-lid top domain boundary and prescribes a fixed scalar concentration
$c_{\textit{top}}$
. On the solid surfaces of the spheres forming the porous medium, a fixed uniform scalar concentration
$c_{\textit{pm}}$
is prescribed. This Dirichlet-type isoscalar boundary condition allows scalar fluxes to pass through the sphere surfaces. As sketched in figure 1, the concentration difference
$\Delta c = c_{\textit{top}} - c_{\textit{pm}}$
generates a downward-directed scalar flux
$\langle \overline {J}_z \rangle$
, which crosses the free-flow region and is then absorbed by the sphere surfaces. This absorption process distinguishes the present study from our previous one (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
), where a strictly zero-flux Neumann boundary condition was applied on the sphere surfaces.
Profiles of the in-plane porosity for different domain extents. The base area
$A_0$
is specified in terms of the sphere diameter
$D$
. The vertical coordinate
$z$
is defined such that the interface
$z=0$
is found where
$\partial ^2 \theta \kern-0.5pt/\kern-1pt \partial z^2 = 0$
. The region
$z/{\kern-1pt}D \gt 2$
is cut off, as
$\theta = 1$
holds on the complete free-flow region.
3.2. Representation of the porous medium
To conduct the flow simulations, mono-disperse random sphere packs with different base areas were generated, as detailed by v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024). Despite their surface irregularity on grain scale, these sphere packs exhibit a level mean bed surface and do not include macroscopic bedforms such as riffles or dunes. The absence of repeating large-scale patterns was confirmed by an evaluation of the spatial autocorrelation function of the bed elevation, which rapidly decays to zero (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025b ).
To conduct simulations with different ratios of flow depth
$h$
and sphere diameter
$D$
, sediment beds with different base area extents in terms of
$D$
are synthesised. For each base area extent, figure 2 shows the in-plane porosity profile
$\theta (z)$
for the sphere pack that is later used in the flow simulation. In analogy to Voermans, Ghisalberti & Ivey (Reference Voermans, Ghisalberti and Ivey2017), the bed-normal coordinate
$z$
is defined such that the interface
$z=0$
coincides with the inflection point
$\partial ^2 \theta \kern-1pt/ \partial z^2 = 0$
of the porosity profile. This purely geometrically defined interface is known before the start of the flow simulation and is therefore used to specify the nominal flow depth
$h$
. Below the strong porosity changes in the near-interface region, a nearly constant bulk porosity of
$\theta _{por} \approx 0.385$
is observed in deeper regions of the sediment bed. Based on the geometric interface, the sediment bed has a depth of
$5 D$
for all simulated cases. The domain sizes are at least
$L_x = 12.8h$
in the streamwise and
$L_y = 6.4h$
in the spanwise direction (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024). The statistics presented in this study were gathered over a minimum of
$T_s u_b \kern-0.5pt/\kern-1pt L_x \gt 10$
, where
$T_s$
is the sampling time period. While this could be shown to ensure statistical convergence for hydraulically rough cases (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
), sampling times of up to
$T_s u_b \kern-0.5pt/\kern-1pt L_x = 33$
were used for transitionally rough- and smooth-wall cases.
3.3. Parameter space
The parameter space explored in the present study is framed using the friction Reynolds number
${\textit{Re}}_\tau$
and the permeability Reynolds number
${\textit{Re}}_{{\kern-1pt}K}$
. The friction Reynolds number, defined as
${\textit{Re}}_\tau = u_\tau h \kern-0.5pt/\kern-1pt \nu$
, mainly characterises the unconfined free-flow further above the sediment bed surface, where we observed similarities between cases with similar
${\textit{Re}}_\tau$
(v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024). Due to computational constraints, the maximum Reynolds number considered is
${\textit{Re}}_\tau = 500$
. The permeability Reynolds number
${\textit{Re}}_{{\kern-1pt}K} = u_\tau \sqrt {K} \kern-0.5pt/\kern-1pt \nu$
incorporates the square root of the permeability
$K$
as an effective length scale for the pore space (Breugem, Boersma & Uittenbogaard Reference Breugem, Boersma and Uittenbogaard2006). By relating this length scale to the viscous length scale
$\delta _\tau = \nu \kern-0.5pt/\kern-1pt u_\tau$
, the permeability Reynolds number
${\textit{Re}}_{{\kern-1pt}K}$
indicates whether or not small-scale turbulent motion can effectively penetrate the pore space. For this reason,
${\textit{Re}}_{{\kern-1pt}K} \ll 1$
represents effectively impermeable conditions and
${\textit{Re}}_{{\kern-1pt}K} \gg 1$
corresponds to the highly permeable regime. According to Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017), values around
${\textit{Re}}_{{\kern-1pt}K} = O(1)$
represent the transitional regime between both extremes. The value of
$K$
was determined via the Darcy relation
$K = - {\mu u_D}{({\partial p}/{\partial x})^{-1}}$
, where
$u_D$
is the Darcy velocity, and agrees well with the prediction by the Kozeny–Carman equation (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024). The simulated cases are distinguished by different values of
$h/D$
, which in turn facilitates the investigation of varying cross-combinations of
${\textit{Re}}_\tau \in [150, 500]$
and
${\textit{Re}}_{{\kern-1pt}K} \in [0.4, 2.8]$
, as visualised in figure 3. In addition to the discussed Reynolds numbers, table 1 lists other nominal parameters of the simulated flow cases. The term nominal parameters is used to emphasise that the parameters refer to the a priori geometrically defined flow depth
$h$
and the corresponding friction velocity
$u_\tau = \sqrt {g_x h}$
.
Overview of the dimensionless parameter space, including reference points from literature. The grey dashed lines represent fixed ratios between the flow depth
$h$
(i.e. boundary layer thickness) and the sphere diameter
$D$
. As reference points, we refer to Breugem et al. (Reference Breugem, Boersma and Uittenbogaard2006), Voermans et al. (Reference Voermans, Ghisalberti and Ivey2017), Shen, Yuan & Phanikumar (Reference Shen, Yuan and Phanikumar2020), and Karra et al. (Reference Karra, Apte, He and Scheibe2023). Figure adapted from v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024).
Overview of nominal case parameters. The variable
$h$
represents the flow depth above the geometrically defined interface,
$D$
is the sphere diameter,
$\Delta x_{i,{\textit{min}}}^+$
describes the side length of the smallest cubic cells near the interface and
$\Delta x_{i,{\textit{max}}}^+$
the side length of the largest cells in the free-flow region. The friction and permeability Reynolds numbers are defined as
${\textit{Re}}_\tau = u_\tau h \kern-0.5pt/\kern-1pt \nu$
,
${\textit{Re}}_{{\kern-1pt}K} = u_\tau \sqrt {K} \kern-0.5pt/\kern-0.5pt \nu$
, respectively, where
$u_\tau = \sqrt {g_x h}$
is the friction velocity,
$K$
the permeability and
$\nu$
the kinematic viscosity. The equivalent sand roughness
$k_s$
is determined as described by v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024).
3.4. Interface position
As the geometrically defined interface is derived a priori from the porosity profile, it lacks a flow-dynamical justification. Thom (Reference Thom1971) and Jackson (Reference Jackson1981) suggested that flow-dynamically more meaningful interface may be determined a posteriori from the drag distribution. v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024) describe an approach to obtain a drag-based interface definition in detail, which we summarise briefly.
We obtain the double-averaged viscous and pressure drag distribution from the vertical gradient of the double-averaged momentum fluxes and the volume source
$g_x$
. Details are provided in Appendix A. The resulting distribution of the superficially double-averaged drag shows a characteristic near-interface peak, which can be associated with the absorption of momentum from the free-flow region. To characterise this free-flow momentum absorption peak by a few parameters and to distinguish it from the Darcy drag, we applied a curve fitting approach with a fitting function
$f ( z, \mu _z, \sigma _z )$
. The parameter
$\mu _z$
represents the centroid of the free-flow momentum absorption peak and we interpret it as the mean interface position, which agrees conceptionally with the roughness centroid of Hantsis & Piomelli (Reference Hantsis and Piomelli2022). Beyond that, the spread
$\sigma _z$
of the drag peak is interpreted as a length scale associated with the interface region. Based on the interface position
$\mu _z$
, a consistent interface-adapted flow depth
$h^\mu = h-\mu _z$
and consistent friction velocity
$u_\tau ^\mu = \sqrt {g_x h^\mu }$
can be defined. Consequently, also interface-adapted Reynolds numbers are obtained (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024), which are slightly smaller than the nominal Reynolds numbers given in table 1.
The different interface definitions open two ways to specify vertical positions. One of them is as
$z/{\kern-1pt}D$
or
$z\kern-0.5pt/{\kern-0.5pt}h$
, where
$z$
refers to the geometrically defined interface. Alternatively, we can use the near-interface coordinate
$\gamma = (z - \mu _z) \kern-0.5pt/\kern-1pt \sigma _z$
or the free-flow coordinate
$\zeta = (z - \mu _z) \kern-0.5pt/\kern-1pt (h - \mu _z)$
, which reflect the drag-based interface position. Both in the near-interface and the free-flow region, we could unveil similarities in the flow field using the interface-adapted quantities for normalisation, while resorting to the coordinates
$\gamma$
and
$\zeta$
(v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024).
3.5. Numerical methods
The applied pore-resolved single-domain DNS avoids model assumptions. In return, it demands that all temporal and spatial scales in both the free-flow region and the pore space are resolved. The massively MPI-parallel simulations were run by our in-house code MGLET (Manhart et al. Reference Manhart, Tremblay, Friedrich, Jenssen, Andersson, Ecer, Satofuka, Kvamsdal, Pettersen, Jenssen, Andersson, Ecer, Satofuka, Kvamsdal, Pettersen, Periaux and Fox2001; Manhart Reference Manhart2004; Peller et al. Reference Peller, Duc, Tremblay and Manhart2006; Sakai et al. Reference Sakai, Mendez, Strandenes, Ohlerich, Pasichnyk, Allalen and Manhart2019), which is optimised for large-scale DNS and large eddy simulation (LES). The incompressible Navier–Stokes equations are solved using an energy-conserving central second-order finite-volume method. For the time integration, MGLET uses an explicit third-order low-storage Runge–Kutta method (Williamson Reference Williamson1980). A multi-grid approach allows a zonal grid refinement, which was applied to ensure that the smallest scales, i.e. Kolmogorov and Batchelor scales, are resolved all across the domain. The near-interface region has the highest demands concerning the spatial resolution. Due to the sphere surfaces, no preferential orientation of boundary layers nor turbulent motion is identifiable, such that we use cubic cells with a side length of
$\Delta x_i^+ \lesssim 1$
(see table 1), which also improves the resolution of the pore geometry. Grid studies showed that a resolution of 48 cells per sphere diameter yields a good resolution of the pore geometries, such that we never go below this resolution anywhere within the sphere pack (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024). MGLET employs an immersed boundary method (Peller et al. Reference Peller, Duc, Tremblay and Manhart2006; Peller Reference Peller2010) to represent the geometry of the complex geometry of the sphere pack. A ghost-cell approach enforces the no-slip boundary condition on the sphere surfaces with second-order spatial accuracy, while the conservation of mass is ensured (Peller Reference Peller2010). The fixed-value Dirichlet boundary condition for the scalar is also implemented with second-order spatial accuracy at the immersed boundary. For a Schmidt number of
${\textit{Sc}} = 1$
, the diffusion boundary thickness is similar to the thickness of the viscous boundary layer, such that the same resolution requirements apply for the scalar and velocity field on the sphere surfaces. Apart from fluxes across the boundaries, strict conservation of the scalar quantity is ensured.
4. Results
The scalar field is analysed (§ 4.1) before the scalar flux and the contributions of different transport processes are addressed (§ 4.2). We describe the absorption of the scalar on the surface of the porous medium (§ 4.3) and make use of it to unveil similarities (§ 4.4). Wherever appropriate, differences and similarities to the findings reported by v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024, Reference v. Wenczowski and Manhart2025a ) are mentioned.
4.1. Scalar concentration field
Figure 4 shows an exemplary realisation of the instantaneous scalar field as well as the time-averaged scalar field within an arbitrarily chosen slice through case M-500. In the free-flow region above the sediment bed, the scalar field appears to be well mixed by eddies of the turbulent fluid motion, as visualised in figure 4(a). The time-averaged field in figure 4(b) shows that a thin layer with high scalar concentration gradient exists near the free-slip top boundary of the domain, where the boundary condition attenuates the vertical fluid motion. In the pore space of the porous medium, the scalar concentration is nearly homogeneous and has a value close to
$c_{\textit{pm}}$
, which is prescribed on the sphere surfaces. The observed spatial homogeneity is in strong contrast to the high spatial variance of the sub-interface scalar field that we observed in a configuration with adiabatic conditions on the sphere surfaces (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
). The dashed boxes in figure 4(b) mark a few individual positions, where larger pores are found near the interface. At these positions, down-welling motion in the time-averaged flow field transports fluid with a higher scalar concentration from the free-flow region into the space between the spheres.
Scalar field within an arbitrarily chosen
$x$
–
$z$
-slice through simulation case M-500. Both the (a) instantaneous and (b) time-averaged scalar fields are normalised by the difference
$\Delta c = c_{\textit{top}} - c_{\textit{pm}}$
, such that the values lie within
$[0,1]$
. Vertical and streamwise coordinates are given in
$x/D$
and
$z/{\kern-1pt}D$
, respectively, where
$D$
is the sphere diameter.
The above-mentioned observations are reflected by the double-averaged scalar concentration profiles in figure 5. In combination with the boundedness of the scalar field, the applied normalisation by
$\Delta c = c_{\textit{top}} - c_{\textit{pm}}$
ensures that the profile spans from zero to unity. In figure 5(a), the bed-normal coordinate
$z$
is normalised by the nominal flow depth
$h$
. A small gradient
$\partial \langle \overline {c} \rangle \kern-0.5pt/\kern-1pt \partial z$
in the centre of the free-flow region hints at a high effective diffusivity with respect to the vertical scalar flux, which is contrasted by larger gradients near the top-boundary of the domain. Near the sediment–water interface around
$z/h=0$
, the scalar concentration profile has a small gradient for cases with high
${\textit{Re}}_{{\kern-1pt}K}$
. With decreasing
${\textit{Re}}_{{\kern-1pt}K}$
, the near-interface gradient increases, such that the profiles gain similarity with the profiles for flow over smooth and impermeable walls. Figure 5(b) uses the sphere diameter to normalise the vertical coordinate
$z$
and focuses on the near-interface region. Below the interface, i.e. at
$z/\!D \lt 0$
, the profiles begin to group according to
${\textit{Re}}_{{\kern-1pt}K}$
, as highlighted by the detail view. At higher
${\textit{Re}}_{{\kern-1pt}K}$
, the normalised scalar concentration profile decays slower with increasing depth. A possible explanation is that high scalar-concentration fluid from the free-flow region is more likely to entrain into the pore space.This motivates an evaluation of the scalar fluxes and the responsible transport processes in the next sections.
Double-averaged scalar profiles plotted against (a)
$z\kern-0.5pt/{\kern-0.5pt}h$
and (b)
$z/{\kern-1pt}D$
, where
$h$
is the nominal flow depth and
$D$
the sphere diameter. The profiles are normalised by the difference
$\Delta c = c_{\textit{top}} - c_{\textit{pm}}$
, such that the values lie within
$[0,1]$
.
4.2. Scalar flux and transport processes
Figure 6(a) uses the sphere diameter
$D$
as a length scale to normalise both the vertical coordinate and the superficially double-averaged scalar flux. It can be seen that the downward-oriented scalar flux remains constant across the free-flow region (
$z/\!D \gt 1$
), whereas the flux starts to decline in amplitude below the crests of the topmost spheres, where scalar absorption occurs. In deeper regions of the sphere pack below
$z/\!D \approx -2$
, hardly any vertical scalar flux is observed, which implies that the incoming flux from the free-flow region has been completely absorbed. Under the given normalisation of the flux with
$\varGamma _c { \Delta c } \kern-0.5pt/\kern-1pt { D }$
, the curves representing cases with similar
${\textit{Re}}_{{\kern-1pt}K}$
collapse in the region shortly below the interface (
$0 \gt z/\!D \gt -2$
), which suggests that the near-interface scalar absorption happens on a length scale comparable to the sphere diameter, while it is controlled by the permeability Reynolds number. In § 4.3, more details on the scalar absorption will follow. Figure 6(b) uses the flow depth
$h$
as the length scale for normalisation. When the scalar flux is normalised with
$\varGamma _c { \Delta c } \kern-0.5pt/\kern-0.5pt { h }$
and plotted against
$z\kern-0.5pt/{\kern-0.5pt}h$
, curves representing cases with similar
${\textit{Re}}_\tau$
group loosely in the free-flow region, whereas no obvious similarity is seen near the interface. The grouping according to
${\textit{Re}}_\tau$
suggests that the friction Reynolds number primarily controls the dimensionless scalar flux across the length scale
$h$
of the free-flow region, whereas a secondary effect seems to remain, which is addressed in the next paragraph.
Profiles of the superficially double-averaged scalar flux
$\theta \langle \overline {J}_z \rangle = \langle \overline {J}_z \rangle ^s$
plotted against (a)
$z/{\kern-1pt}D$
and (b)
$z\kern-0.5pt/{\kern-0.5pt}h$
, where
$D$
is the sphere diameter and
$h$
is the flow depth. The flux itself is normalised by a diffusive flux, which uses either
$D$
or
$h$
as a length scale, in addition to the scalar diffusion coefficient
$\varGamma _c$
and the scalar concentration difference
$\Delta c = c_{\textit{top}} - c_{\textit{pm}}$
.
In the following, the constant double-averaged scalar flux in the free-flow region will be referred to as
$J_{\textit{ff}}$
. The ratio of
$J_{\textit{ff}}$
and
$\varGamma _c { \Delta c } \kern-0.5pt/\kern-1pt { h }$
can be interpreted as a Sherwood number. The Sherwood number
${\textit{Sh}}$
for scalar mass transfer is equivalent to the Nusselt number
$Nu$
for heat transfer, as it puts the actual scalar flux into relation to the expected flux under purely diffusive conditions. As a primary trend, figure 7 shows that the Sherwood number increases with increasing
${\textit{Re}}_\tau$
. Among cases with similar
${\textit{Re}}_\tau$
, the flux increases slightly with increasing
${\textit{Re}}_{{\kern-1pt}K}$
, as it is best visible for the cases with nominal
${\textit{Re}}_\tau = 300$
. This secondary trend implies that higher roughness and permeability facilitate the scalar transfer from the bed to the top surface of the domain.
Sherwood number
${\textit{Sh}}$
plotted over the nominal friction Reynolds number
${\textit{Re}}_\tau$
. Among cases with similar
${\textit{Re}}_\tau$
, the Sherwood number increases with increasing
${\textit{Re}}_{{\kern-1pt}K}$
.
Equation (2.9) allows a decomposition of the total scalar flux
$ \langle \overline {J}_z \rangle ^s$
within the fluid into contributions of turbulent transport, dispersive transport and diffusive transport. Figure 8(a) shows the total flux under normalisation with
$J_{\textit{ff}}$
. When plotted against
$z/{\kern-1pt}D$
, the curves group according to
${\textit{Re}}_{{\kern-1pt}K}$
. With increasing
${\textit{Re}}_{{\kern-1pt}K}$
, the total scalar flux entrains deeper into the sediment bed, before it is absorbed. At higher
${\textit{Re}}_{{\kern-1pt}K}$
, the wall-blocking effect is relaxed and turbulent motion can penetrate the near-interface pore space (e.g. v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024), which explains the more intense turbulent scalar transport near the interface shown in figure 8(b). Even at
${\textit{Re}}_{{\kern-1pt}K} = 2.8$
, however, the turbulent transport has decreased to a value of approximately
$0.05 \, J_{\textit{ff}}$
at a depth of
$z/\!D = -0.5$
. The profiles in figure 8(c) show that dispersive scalar transport is nearly negligible above the crests of the spheres, while it peaks between
$z/\!D = 0$
and
$z/\!D=0.2$
. The amplitude of the peak grows with increasing
${\textit{Re}}_{{\kern-1pt}K}$
, while the position of the peak lowers. Below the interface, the dispersive transport declines more slowly than the turbulent transport. At
$z/\!D = -0.5$
, the dispersive transport accounts for more than half of the total flux for all cases with
${\textit{Re}}_{{\kern-1pt}K} \gt 1$
. The scalar transport due to molecular diffusion peaks slightly above the interface, as shown by figure 8(d). The peak is most emphasised for the cases with low
${\textit{Re}}_{{\kern-1pt}K}$
, for which we observed steeper gradients in the double-averaged scalar concentration profile (see figure 5). At a depth of
$z/\!D = -0.5$
, the diffusive scalar flux has nearly entirely decayed.
Near-interface profiles of double-averaged scalar fluxes under normalisation with
$J_{\textit{ff}}$
, i.e. the total scalar flux in the free-flow region. The normalised profiles of the (a) total, (b) turbulent, (c) dispersive and (d) diffusive transport group according to
${\textit{Re}}_{{\kern-1pt}K}$
. The vertical coordinate
$z$
is normalised by the sphere diameter
$D$
.
To gain further insight, we analyse the spatial distributions of the three transport processes. Figure 9 shows an arbitrarily chosen slice through case M-500. The local contributions to turbulent, dispersive and diffusive scalar transport are presented under normalisation with
$J_{\textit{ff}}$
. As expected, turbulent transport dominates in the free-flow region. Although not perfectly, the distribution in the outer region is relatively homogeneous around
$ \overline { w^\prime c^\prime } \kern-0.5pt/\kern-1pt J_{\textit{ff}} \approx 1$
. Spatial variations are visible directly above the sphere surfaces and protrude to approximately
$2D$
into the free flow. The dashed boxes mark regions, where large pores near the interface allow turbulent transport, whereas, apart from those locations, hardly any turbulent transport occurs within the pore space. Diffusive scalar transport dominates in a thin layer near the top domain boundary, where the applied boundary conditions suppress advective scalar transport in the vertical direction. Thin diffusion-dominated boundary layers also form at the crests of the topmost spheres, where strong scalar gradients between the isoscalar no-slip sphere surface and the fluid in the well-mixed free-flow region are expected. Dispersive transport is primarily observed near the topmost layer of spheres. Its distribution in space is not homogeneous, which hints at an influence of local features of the sediment bed. In particular, the individual larger near-interface pores within the dashed boxes act as hotspots of dispersive transport. A purely visual analysis suggests that these larger pores tend to be located on the windward side of protruding spheres, where locally higher hydrodynamic pressure is expected. Due to the locally higher pressure, some of the oncoming fluid is pushed upwards, whereas some is also pushed into the pore space. This explains the downwelling fluid motion, which advects fluid with higher scalar concentration into the pore space (see also figure 4).
Local contributions to (a) turbulent, (b) dispersive and (c) diffusive scalar transport within an arbitrarily chosen
$x$
–
$z$
-plane of simulation case M-500. The values are normalised by
$J_{\textit{ff}}$
, i.e. by the total scalar flux in the free-flow region. Coordinates are given in
$x/D$
and
$z/{\kern-1pt}D$
, where
$D$
is the sphere diameter.
The above observation is a qualitative difference to a configuration with zero-flux boundary conditions on the sphere surfaces, where upwelling regions (‘chimneys’) induce spatial variations in the scalar field by transporting fluid with different scalar concentration from deeper regions to the upper layer of the sediment (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a ).
4.3. Scalar absorption
Previously (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024), we documented that the drag-based interface (see § 3.4) can be used as a framework for normalisation and coordinates in which similarities of flow quantity profiles can be found. This motivates a closer look at the scalar absorption. Under steady conditions, (2.8) and (2.10) allow us to compute the superficially double-averaged scalar absorption
$f_a^s$
, i.e. the diffusive scalar flux from the fluid to the solid, from the vertical gradient of the superficial scalar flux inside the fluid domain,
\begin{equation} f_a^s \: = \: {\textit{BT}}_{j}^{\,s}\left ( \overline {J}_{\!j} \right ) \: = \: \frac {1}{A_0} \oint _{s} \, - \, \varGamma _c \frac {\partial \overline {c}}{\partial x_{\!j}} \, \frac {n_{\!j}}{\sqrt {n_x^2 + n_y^2}} \, {\rm d}s \: = \: - \frac {\partial }{\partial z} \left ( \theta \!\left \langle \overline {J}_z \right \rangle \right )\! . \end{equation}
The distributions of the scalar absorption are shown in figure 10(a). The normalisation by
$J_{\textit{ff}}$
implies that the area under the curve must be unity, as the total scalar flux from the free-flow region must be absorbed. For cases with comparable
${\textit{Re}}_{{\kern-1pt}K}$
, the scalar absorption distributions are similar. With increasing
${\textit{Re}}_{{\kern-1pt}K}$
, the amplitude of the peak decreases and the distributions spread and reach to greater depths, while the peak position moves downward. Figure 10(b) provides the distributions of the superficially double-averaged momentum absorption. In comparison to the scalar absorption, the near-interface peaks in the momentum absorption distributions tend to have a higher amplitude and a smaller spread. In contrast to
$f_a^s$
, the momentum absorption at greater depth is not zero, as the volume force
$g_x$
acts as a momentum source on the fluid within the pore space.
Distribution of the superficially double-averaged (a) scalar absorption
$\textstyle f_a^s$
and (b) momentum absorption
$ \textstyle f_{(p+\nu )}^s$
within the near-interface region. For normalisation, the total scalar flux in the free-flow region
$\textstyle J_{\textit{ff}}$
, the nominal friction velocity
$\textstyle u_\tau$
and the sphere diameter
$D$
are used.
To quantify the observed trends, we evaluate statistical moments of the scalar absorption distribution, which are summarised in table 2. The notation with a hat, e.g.
$\hat {\mu }_{z}$
and
$\hat {\sigma }_{z}$
, emphasises that values refer to the scalar absorption distribution, while
$\mu _z$
and
$\sigma _z$
without a hat are based on the drag distribution (see § 3.4). To make trends clearly visible, figure 11 provides the mean value
$\hat {\mu }_{z}$
and the standard deviation
$\hat {\sigma }_{z}$
over
${\textit{Re}}_{{\kern-1pt}K}$
, while
$\mu _{z}$
and
$\sigma _{z}$
are plotted as grey symbols for direct comparison. Figure 11(a) shows that the mean position of scalar absorption,
$\hat {\mu }_{z}$
, lowers with increasing
${\textit{Re}}_{{\kern-1pt}K}$
. A similar, yet less emphasised, behaviour is observed for
$\mu _{z}$
. In general, we observe that
$\hat {\mu }_z \lt \mu _{z}$
, which implies that the scalar is absorbed at a greater depth than the momentum. The standard deviations
$\hat {\sigma }_{z}$
and
$\sigma _{z}$
are shown in figure 11(b) for direct comparison. As
$\hat {\sigma }_{z} \gt \sigma _{z}$
, the scalar absorption spreads over a wider range in terms of
$z$
than the momentum absorption, which appears to happen within a comparatively thin layer. Figure 11(c) shows the skewness of the scalar absorption distribution over
${\textit{Re}}_{{\kern-1pt}K}$
. The amplitude of the negative skewness increases with progressively higher
${\textit{Re}}_{{\kern-1pt}K}$
. This trend, however, seems to saturate for
${\textit{Re}}_{{\kern-1pt}K} \gt 2$
. The skewness parameter is not available for the momentum absorption. Together, the comparisons suggest that momentum is absorbed more effectively by the porous medium than the scalar, which leads to recognisable differences between the absorption distributions. Higher
${\textit{Re}}_{{\kern-1pt}K}$
seem to emphasise those differences.
Description of the near-interface scalar absorption distribution by statistical moments of different order. The values are normalised by the sphere diameter
$D$
.
Parameters describing the scalar absorption distribution, i.e. (a) mean position
$\hat {\mu }_z$
, (b) spread, quantified by the standard deviation
$\hat {\sigma }_z$
, and (c) skewness
$skew$
of the distribution. In panels (a) and (b), grey symbols represent the corresponding parameters of the momentum absorption distribution (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024). Where appropriate, values are normalised by the sphere diameter
$D$
.
4.4. Near-interface similarities
Previously (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024), we used the mean position and the spread of the momentum absorption to define an interface coordinate, which exposed similarities in the near-interface velocity field. That motivates us to use the parameters
$\hat {\mu }_z$
and
$\hat {\sigma }_z$
of the scalar absorption distribution similarly to define a scalar near-interface coordinate
$\hat {\gamma }$
as
\begin{equation} \hat {\gamma } = \frac { ( z - \hat {\mu }_z ) }{ \hat {\sigma }_z } . \end{equation}
In addition to the interface coordinate in (4.2), we define the friction scalar concentration as
$c_\tau = { J_{\textit{ff}} } \, \kern-0.5pt/\kern-1pt \, { u_\tau ^\mu }$
, which in turn comprises the friction velocity
$u_\tau ^\mu = \sqrt { g_x (h - \mu _z) }$
.
In contrast to figure 12(a), which uses
$z/{\kern-1pt}D$
on the vertical axis, figure 12(b) shows the profiles of the normalised scalar flux plotted against
$\hat {\gamma }$
. The curves collapse with good accuracy, which confirms that the case-specific parameters
$\hat {\mu }_z$
and
$\hat {\sigma }_z$
alone describe the scalar absorption distribution sufficiently well. Only a detailed view with strong zoom-in on the region
$\hat {\gamma } \in [-3,-1]$
shows small differences in the tailing behaviour of the profiles, which can be associated with the differences in the skewness or other higher order moments among the scalar absorption distributions. In figures 12(c) and 12(d), the double-averaged scalar concentration profiles normalised by the friction scalar concentration
$c_\tau$
are plotted against
$z/{\kern-1pt}D$
and
$\hat {\gamma }$
, respectively. If the interface coordinate
$\hat {\gamma }$
is used, a similarity can be seen, whereas solely case S-150 deviates slightly from the remaining profiles, as emphasised by the embedded detail view. The transitionally rough regime of case S-150 may explain the differences to the other cases, which can be categorised as hydraulically rough. In contrast to the collapse of the flux profiles, the partial collapse of the scalar profiles is not directly connected to the definition of the absorption-based coordinate
$\hat {\gamma }$
, as it also includes the effect of near-interface scalar transport.
Near-interface similarities of the normalised double-averaged, scalar flux and scalar concentration. In panels (a) and (c), the quantities are plotted against
$z/{\kern-1pt}D$
, where
$z=0$
is defined by the inflection point of the porosity profile and
$D$
is the sphere diameter. In panels (b) and (d), the interface coordinate
$\hat {\gamma }$
is used, which considers the mean position
$\hat {\mu }_z$
and the spread
$\hat {\sigma }_z$
of the scalar absorption. The scalar flux is normalised by the flux
$J_{\textit{ff}}$
in the free-flow region. The friction scalar concentration
$c_\tau$
is used to normalise the scalar concentration profile. The concentration
$c_{\textit{pm}}$
is prescribed on the porous medium.
The above-mentioned observations encourage us to analyse if the temporal scalar fluctuations
$\langle \overline {c' c'} \rangle$
and the spatial scalar variations
$\langle \widetilde {\overline {c}} \, \widetilde {\overline {c}} \rangle$
also exhibit similarities within the interface coordinate framework. By plotting the normalised scalar fluctuations against
$z/{\kern-1pt}D$
, figure 13(a) shows that scalar fluctuations can be observed at progressively greater depth below the geometrically defined interface as
${\textit{Re}}_{{\kern-1pt}K}$
increases. Figure 13(b) refers to the interface coordinate
$\hat {\gamma }$
and shows similarities among the profiles of all hydraulically rough cases, whereas again, the transitionally rough case S-150 renders an exception. Like in figure 12(d), the profile of case S-150 is shifted upwards in comparison to the collapsing profiles of the remaining cases. In figure 13(c), the spatial variations in the time-averaged scalar field are plotted against
$z/{\kern-1pt}D$
. The spatial variations peak near the interface and penetrate to progressively greater depth below the geometrically defined interface with increasing
${\textit{Re}}_{{\kern-1pt}K}$
. When plotted against
$\hat {\gamma }$
, the positions of the scalar variation peaks collapse for all hydraulically rough cases, as shown in figure 13(d). The profile for the transitionally rough case S-150 shows a qualitatively different behaviour, as both peak amplitude and position deviate considerably from the other cases.
Near-interface similarities of the normalised temporal scalar fluctuations and spatial scalar variations. In panels (a) and (c), the quantities are plotted against
$z/{\kern-1pt}D$
, where
$z=0$
is defined by the inflection point of the porosity profile and
$D$
is the sphere diameter. In panels (b) and (d), the interface coordinate
$\hat {\gamma }$
is used, which considers the mean position
$\hat {\mu }_z$
and the spread
$\hat {\sigma }_z$
of the scalar absorption. The friction scalar concentration
$c_\tau$
is used for normalisation.
In summary, the scalar absorption distribution appears to be helpful when it comes to finding near-interface similarities of scalar-related quantities. As such, it plays a similar role that the drag distribution plays for similarities in the flow field.
4.5. Outer-layer similarity
While the previous sections focused on the interface region, we address the region above the roughness sublayer in the following. Previously (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
), we used the defect representation to show similarities in the dimensionless scalar profiles of cases with comparable
${\textit{Re}}_\tau$
in the region far from the sediment bed surface. Those similarities were linked to the observation that the far-wall effective diffusivity depends on
${\textit{Re}}_\tau$
. Figure 14 provides the dimensionless scalar profiles
$c^+ = ( \langle \overline {c} \rangle - c_{\textit{pm}} ) \, \kern-0.5pt/\kern-1pt \, c_\tau$
of cases with comparable
${\textit{Re}}_\tau$
. Note that the profiles are plotted over a dimensionless coordinate, which has zero value at the drag-based interface
$z=\mu _z$
and a unity value at the top domain boundary. For each group of cases with comparable
${\textit{Re}}_\tau$
, a smooth-wall case serves as a reference, such that the scalar roughness function
$\Delta c^+$
for the cases with
${\textit{Re}}_{{\kern-1pt}K}\gt 0$
can be computed. By down-shifting the smooth-wall reference profile by
$\Delta c^+$
, the dotted lines allow a check of the profile similarity.
Dimensionless scalar profiles over the complete flow depth under normalisation with
$c_\tau$
. The horizontal axis represents the bed-normal direction, where
$\mu _z$
marks the drag-based interface position. For each group of cases with similar
${\textit{Re}}_\tau$
, the smooth-wall case acts as a reference for the computation of
$\Delta c^+$
-values, which are given in braces for each case. The dotted lines result from a downward shift of the smooth-wall profile by
$\Delta c^+$
and show the similarity in the far-wall region.
The scalar roughness function deserves some further comments. Figure 15(a) shows that
$\Delta c^+$
increases monotonously over
${\textit{Re}}_{{\kern-1pt}K}$
. Figure 15(b) shows the relation to the roughness function values
$\Delta u^+$
, which were determined by v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024) and describe the deviation from a smooth-wall velocity profile at comparable
${\textit{Re}}_\tau$
. Both roughness functions behave similarly, such that the data points lie near the reference line defined by
$\Delta c^+ = \Delta u^+$
, whereas
$\Delta u^+$
is slightly larger than
$\Delta c^+$
. For impermeable rough-wall cases at
${\textit{Sc}} = 1$
, the data gathered by Hantsis & Piomelli (Reference Hantsis and Piomelli2024) suggest a more emphasised tendency towards
$\Delta c^+ \lt \Delta u^+$
.
Behaviour of the scalar roughness function
$\Delta c^+$
including (a) dependency on the permeability Reynolds number
${\textit{Re}}_{{\kern-1pt}K}$
and (b) relation to the velocity roughness function
$\Delta u^+$
. The dotted line marks
$\Delta c^+ = \Delta u^+$
as a reference.
The above-mentioned observations motivate an evaluation of the Reynolds analogy factor
${\textit{RA}}$
. Following MacDonald et al. (Reference MacDonald, Hutchins and Chung2019) and Hantsis & Piomelli (Reference Hantsis and Piomelli2024), we define the Reynolds analogy factor as
\begin{equation} \textit{RA} = \frac {u_b^+}{c_m^+} \quad \text{with} \quad u_b = \frac {1}{h-\mu _z} \int _{\mu _z}^h \langle \overline {u} \rangle \, {\rm d}z \; , \; c_m = \frac { \int _{\mu _z}^h (\langle \overline {c} \rangle - c_{\textit{pm}}) \langle \overline {u} \rangle \, {\rm d}z \: }{ \int _{\mu _z}^h \langle \overline {u} \rangle \, {\rm d}z } \; . \end{equation}
In (4.3),
$u_b$
represents the bulk velocity and
$c_m$
the mixed-mean scalar concentration. Through the definition of
$c_m$
, the Reynolds analogy factor not only depends on the shifts
$\Delta u^+$
and
$\Delta c^+$
, but also accounts for the shape of the velocity and scalar profiles. Accordingly, the parameter is also influenced by
${\textit{Re}}_{{\kern-1pt}K}$
-dependent changes of the profiles in the near-interface region, where the profiles deviate from a mere shift of the smooth-wall profile. Figure 16 shows that the normalised Reynolds analogy factor declines over
${\textit{Re}}_{{\kern-1pt}K}$
, which agrees well with the trend reported for rough-wall flow (Hantsis & Piomelli Reference Hantsis and Piomelli2024). For the normalisation, the Reynolds analogy factor of the smooth-wall simulation at comparable
${\textit{Re}}_\tau$
is used. Note that the reference values are less than unity as
${u_b^+} \lt {c_m^+}$
, which is a result of the different boundary conditions for scalar and velocity at the top of the domain as well as of the absence of a volume source for the scalar. Due to this initial inequality
${u_b^+} \lt {c_m^+}$
, even a decrease of both quantities by the same amount leads to a decrease of the (normalised) Reynolds analogy factor, rendering the parameter less meaningful under the present boundary conditions.
Normalised Reynolds analogy factors
${\textit{RA}}$
as a function of
${\textit{Re}}_{{\kern-1pt}K}$
. The reference value
${\textit{RA}}_{\textit{ref}}$
for the normalisation corresponds to the
${\textit{RA}}$
of the smooth-wall case at comparable
${\textit{Re}}_\tau$
.
Budget of scalar and velocity fluctuations within the near-interface region. The terms in (a) the budget of
$\langle \overline {c' c'} \rangle /2$
and the terms in (b) the budget of
$\langle \overline {u_{\!j}^{\prime} u_{\!j}^{\prime}} \rangle /2$
are normalised by the sphere diameter
$D$
, the friction velocity
$u_\tau ^\mu$
and the friction scalar concentration
$c_\tau$
. Case M-500 was used as an example.
4.6. Similarity of scalar variance and TKE budgets
In the following, we focus on the near-interface production mechanisms of scalar fluctuations and compare them with the processes involved in the production of velocity fluctuations, i.e. production of turbulent kinetic energy (TKE). Using case M-500 with
${\textit{Re}}_{{\kern-1pt}K} = 2.8$
as an example, figure 17(a) shows the budget for
$\langle \overline {c' c'} \rangle /2$
under normalisation by
$u_\tau ^\mu$
,
$c_\tau$
and
$D$
. The total production
$\varPi _{{tot}}$
reaches its peak value close to the sphere crests at
$z/\!D \approx 0.5$
. The gradient production
$\varPi _{\textit{grad}}$
seems to be the main process responsible for this peak. The peak value of the form-induced production is found near
$z/\!D \approx 0$
and its amplitude is critically smaller than the maximal gradient production value. Only at depths below
$z/\!D \lesssim 0.5$
, the form-induced production accounts for the larger part of the production of scalar fluctuations. The total transport
$T_{\textit{tot}}$
relocates scalar fluctuations from the position of the maximal
$\varPi _{{tot}}$
further downwards into the pore space. Figure 17(b) shows the TKE budget for the same case under normalisation with
$u_\tau ^\mu$
and
$D$
. The direct comparison shows that the production maxima of scalar fluctuations and TKE share a similar amplitude and are found approximately at the same position. Mainly the shear production, which corresponds to the gradient production of
$\langle \overline {c' c'} \rangle /2$
, contributes to the peak, whereas form-induced production only dominates at greater depth. Compared with the scalar fluctuations budget, the downward-oriented transport of TKE is slightly more intense due to pressure transport (e.g. v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025b
). The parallels between the near-interface budgets of scalar fluctuations and TKE, particularly in the production terms, motivate a further investigation of the correlation between scalar and velocity fluctuations.
Plane-averaged correlation coefficients between scalar fluctuations
$c'$
and (a) streamwise velocity fluctuations
$u'$
and (b) bed-normal velocity fluctuation
$w'$
within the near-interface region. The curves group according to
${\textit{Re}}_{{\kern-1pt}K}$
.
4.7. Scalar-velocity correlations
The previous section has shown parallels between the production mechanisms of scalar fluctuations and velocity fluctuations. As the Reynolds analogy implies that scalar and momentum fluxes are driven by the same transport mechanisms, the turbulent fluxes
$\overline {u'w'}$
and
$\overline {c'w'}$
should also exhibit certain parallels, which may weaken with
${\textit{Re}}_{{\kern-1pt}K}$
as suggested by figure 16. Those observations motivate an investigation of the correlation between
$c'$
and the velocity fluctuations
$u'$
as well as
$w'$
. Accordingly, the plane-averaged correlation coefficients are defined as
\begin{equation} \left \langle C_{c'u'} \right \rangle = \left \langle \frac {\overline {c'u'}}{\sqrt {\overline {c'c'}} \, \sqrt {\overline {u'u'}}} \right \rangle \hspace {1cm} \text{and} \hspace {1cm} \left \langle C_{c'w'} \right \rangle = \left \langle \frac {\overline {c'w'}}{\sqrt {\overline {c'c'}} \, \sqrt {\overline {w'w'}}} \right \rangle \!. \end{equation}
Both correlation coefficients in (4.4) can have values
$\in [-1,1]$
, where positive values indicate a positive linear correlation and negative values a negative linear correlation.
Figure 18(a) shows a clear positive correlation between scalar and streamwise velocity fluctuations for
$z/\!D \gt -0.2$
. For the high-
${\textit{Re}}_{{\kern-1pt}K}$
cases L-300 and M-500, the correlation coefficient reaches its maximum at a value of
$\langle C_{c'u'} \rangle \approx 0.6$
. In contrast to that, the strongest linear correlation with
$\langle C_{c'u'} \rangle \approx 0.9$
is observed for the transitionally rough case S-150. Similar values are reported by Pirozzoli et al. (Reference Pirozzoli, Bernardini and Orlandi2016) for smooth-wall cases, where the strongest correlation is observed in the buffer layer. Below
$z/\!D \approx -0.2$
, the correlation has negative values of small amplitudes, which may be a result of counter-streamwise scalar transport due to recirculating flow in this region. In the complete near-interface region, the curves group according to
${\textit{Re}}_{{\kern-1pt}K}$
. The correlation coefficient
$\langle C_{c'w'} \rangle$
is shown in figure 18(b). Negative values of the coefficient can be associated with the downward-oriented turbulent scalar transport. Near the interface and below, the curves group according to
${\textit{Re}}_{{\kern-1pt}K}$
, whilst their amplitudes decline slower for cases with higher
${\textit{Re}}_{{\kern-1pt}K}$
. In the region
$z/\!D \gt 0.5$
, all cases show values of
$\langle C_{c'w'} \rangle \approx -0.4$
such that the correlation between
$w'$
and
$c'$
in this region is weaker than the strong correlation between
$u'$
and
$c'$
.
By showing realisations of
$c'/c_\tau$
and
$u'/u_\tau ^\mu$
at the same time side-by-side, figure 19 provides a visual impression of the similarity between scalar and streamwise velocity fluctuations. For the transitionally rough case S-150, the maximum of
$\langle C_{c'u'} \rangle$
is found at
$z/\!D = 0.79$
, where a streaky pattern can be seen both in the streamwise velocity fluctuation and in the scalar fluctuation. The strong correlation and the sharp contours facilitate the identification of corresponding structures in
$c'$
and
$u'$
. In the fully rough case M-500,
$\langle C_{c'u'} \rangle$
peaks at
$z/\!D = 0.56$
. Whereas steep gradients of
$c'$
lead to strong colour contrasts in the scalar fluctuation field, the field
$u'$
appears blurry and poorer in contrast, which is also due to lower amplitudes of extreme values. The weaker correlation makes it harder to detect a one-to-one correspondence among the individual structures in the velocity and scalar fields. With increasing
${\textit{Re}}_{{\kern-1pt}K}$
, streamwise streaks in both
$c'$
and
$u'$
seem to vanish, such that bulkier structures dominate the pattern of both fields.
Fluctuations of scalar and streamwise velocity field in comparison. For case (a) S-150 and (b) M-500, the normalised realisations of
$c'$
and
$u'$
were captured at
$z/\!D=0.79$
and
$z/\!D=0.56$
, respectively, which are the positions of the maximal correlation coefficient
$\langle C_{c'u'} \rangle$
(see figure 18). Coordinates are given as
$x^+ = x u_\tau ^\mu \kern-0.5pt/\kern-1pt \nu$
and
$y^+ = y u_\tau ^\mu \kern-0.5pt/\kern-1pt \nu$
. For each case, the images show the same patch at the same time.
5. Discussion
Even at a unitary Schmidt number of
${\textit{Sc}} = \nu \kern-0.5pt/\kern-1pt \varGamma _c = 1$
, a direct comparison shows that the mean position of scalar absorption
$\hat {\mu }_z$
is found at a lower position than the mean position of momentum absorption
$\mu _z$
(see v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024). Also, the scalar absorption spreads over a wider vertical region than the momentum absorption (figure 10). For flow over impermeable rough walls, a widely accepted explanation (e.g. Chung et al. Reference Chung, Hutchins, Schultz and Flack2021; Hantsis & Piomelli Reference Hantsis and Piomelli2024) is that the momentum transfer between the fluid and the solid spheres happens via both viscous and pressure drag. As viscous drag can be interpreted as a diffusion of momentum across the solid–fluid interface, it finds its analogy in the diffusive scalar transport across the interface. An analogy to the pressure drag does not exist for the scalar absorption, such that the pressure drag represents an additional pathway for the momentum transfer, which breaks the Reynolds analogy (Reynolds Reference Reynolds1961). The role of pressure drag (see also Appendix A) would explain the observation that differences between momentum and scalar absorption become more emphasised with increasing
${\textit{Re}}_{{\kern-1pt}K}$
(figure 11) and, thus, also with increasing roughness Reynolds number
$k_s^+$
(table 1). Beyond those parallels to rough-wall flow, increasingly permeable wall behaviour at
${\textit{Re}}_{{\kern-1pt}K} \gt 1$
causes a progressively longer tail of the scalar absorption distribution, which extends into greater depths of the sphere pack and, thus, contributes to the negative skewness of the scalar absorption distribution (figure 11). A decomposition of the scalar flux (figure 8) shows that the scalar is transferred to greater depth mainly by turbulent and dispersive transport due to bed-normal fluid motion, which would be blocked by an impermeable rough wall.
Both scalar transport and momentum transfer in the well-mixed free-flow region are mainly due to turbulent fluid motion, such that Townsend’s wall similarity hypothesis (Townsend Reference Townsend1976) is expected to apply also for the scalar transport. Previously (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2024), we demonstrated that outer-layer similarities in the flow field for cases with comparable
${\textit{Re}}_\tau$
can be seen most clearly under normalisation with the flow depth
$h^\mu$
and the friction velocity
$u_\tau ^\mu$
, i.e. with interface-adapted quantities. Using the interface-adapted friction concentration
$c_\tau = J_{\textit{ff}} /\kern-1pt u_\tau ^\mu$
for normalisation, we could show that the far-wall scalar profile can be replicated via a shift of the smooth-wall profile by a uniform value of
$\Delta c^+$
. For the flow over random sphere packs, we observe that
$\Delta c^+ \approx \Delta u^+$
(figure 15), whereas
$\Delta c^+ \lt \Delta u^+$
is reported for rough-wall cases (e.g. Hantsis & Piomelli Reference Hantsis and Piomelli2024). Larger values of
$\Delta u^+$
and
$\Delta c^+$
can be associated with enhanced absorption of momentum and scalar, respectively. The increase of the share of the pressure drag to the total drag would suggest that the ratio of
$\Delta c^+$
to
$\Delta u^+$
declines with increasing
${\textit{Re}}_{{\kern-1pt}K}$
, which is not the case. A possible explanation for
$\Delta c^+ \approx \Delta u^+$
for our case could lie in the reduced wall-blocking effect with increasing
${\textit{Re}}_{{\kern-1pt}K}$
. As discussed, fluid motion across the interface additionally activates deeper layers for scalar transfer, which increases the active fluid–solid interface area. This effect would enhance scalar transfer and, thus, increase the scalar roughness function
$\Delta c^+$
.
In addition to the turbulent scalar transport in the free-flow region, it has been discussed by several authors that the thickness
$\delta _c$
of the diffusive sublayer influences the dimensionless mass or heat transfer. For the smooth-wall case (
${\textit{Re}}_{{\kern-1pt}K} \ll 1$
), the thickness of the diffusive sublayer has been expressed in the form
$\delta _c \kern-0.5pt/\kern-1pt \delta _\nu \propto Sc^m$
with
$m\in [-0.5, -0.29]$
(Kader Reference Kader1981; Na, Papavassiliou & Hanratty Reference Na, Papavassiliou and Hanratty1999; Schwertfirm & Manhart Reference Schwertfirm and Manhart2007; Alcántara-Á vila & Hoyas Reference Alcántara-Á vila and Hoyas2021). The thickness of the viscous sublayer
$\delta _\nu$
goes with
$\nu /u_\tau$
, and hence,
$\delta _\nu /h\propto Re_\tau ^{-1}$
. Therefore,
$\delta _c \kern-0.5pt/\kern-1pt h \propto Re_\tau ^{-1} Sc^m$
for smooth and impermeable walls. For
${\textit{Re}}_{{\kern-1pt}K} \gg 1$
, however, one can argue that the flow penetrates into the pore space and advects the scalar deeper into the sediment. For this highly permeable limit case, viscous and diffusive boundary layers can develop on the sphere surfaces and the scalar absorption area increases. The exact scaling of those boundary layers will depend on flow type in the sediment – laminar or turbulent – and the Schmidt number. Most likely, this scaling would differ from the smooth-wall scaling. Despite a similarity in the well-mixed free-flow region, the scalar boundary condition on the sphere surfaces strongly influences the scalar field within the porous medium as well as the processes at the sediment bed surface. Previously (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
), we had applied zero-flux Neumann-type scalar boundary conditions on the sphere surfaces, which represents an adiabatic behaviour. Only at the bottom and top domain boundary, fixed scalar concentrations were prescribed to induce a bed-normal scalar flux. Though we also considered
${\textit{Sc}} = 1$
, the scalar field within the adiabatic sphere pack showed strong spatial variations on length scales of several sphere diameters, which became progressively more emphasised with increasing
${\textit{Re}}_{{\kern-1pt}K}$
. This observation is contrasted by the spatial homogeneity of the scalar distribution inside the sphere pack, which even persists at
${\textit{Re}}_{{\kern-1pt}K} \approx 2.8$
under the isothermal boundary condition in the present study (figure 4). The described difference can be explained by different diffusion length scales and, thus, also different time scales for the scalar diffusion process. While in the present study,
$\sqrt {K}$
can be taken as diffusion length scale as the scalar diffuses to the next sphere surface, in the adiabatic case, the scalar variances occur at larger length scales of the order of the bed depth.
Under isothermal sphere surface conditions, the quick diffusion suppresses spatial variations in the scalar field. As a consequence, dispersive transport decays within few sphere diameters below the sediment bed surface (figure 8) and does not affect critically deeper regions, as it was observed under adiabatic conditions (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
). Further, the existence of spatial variations
$\widetilde {\overline {c}}$
in the near-interface region is a prerequisite for the form-induced production of temporal fluctuations in the scalar field, as discussed in § 2.4. That explains why form-induced production plays a major role in the presence of adiabatic scalar boundary conditions (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
), while it plays a secondary role for the isothermal conditions of the present study (figure 17).
Though the evolution of velocity and the scalar field is described by different equations, the values of both fields on the sphere surfaces are prescribed by means of Dirichlet-type boundary conditions. This similarity of the boundary conditions seems to prepare the ground for several parallels between the velocity field and the scalar field. Absorption processes appear to be highly relevant in the near-interface region and can be used to define frameworks, in which near-interface similarities can be seen (e.g. figure 12). Both the double-averaged scalar profile and the streamwise velocity profile exhibit steep gradients in the near-interface region. Therefore, gradient or shear production are the main sources of the near-interface production of scalar fluctuations and TKE, respectively.
The above-mentioned parallels between velocity and scalar field are also reflected by the turbulent Schmidt number
${\textit{Sc}}_t$
, which is a frequently used parameter for modelling purposes, as it represents the ratio between turbulent effective viscosity
$\nu _t = \langle \overline {u' w'} \rangle \kern-0.5pt/\kern-1pt \langle \partial \overline {u} \kern-0.5pt/\kern-1pt \partial z \rangle$
and the turbulent effective diffusivity
$\varGamma _t = \langle \overline {c' w'} \rangle \kern-0.5pt/\kern-1pt \langle \partial \overline {c} \kern-0.5pt/\kern-1pt \partial z \rangle$
. Accordingly, a turbulent Schmidt number for bed-normal scalar transport (MacDonald et al. Reference MacDonald, Hutchins and Chung2019) can be defined as
\begin{equation} Sc_t = \frac {\nu _t}{\varGamma _t} = \frac {\langle \overline {u' w'} \rangle \, \langle \partial \overline {c} \kern-0.5pt/\kern-1pt \partial z \rangle }{\langle \overline {c' w'} \rangle \, \langle \partial \overline {u} \kern-0.5pt/\kern-1pt \partial z \rangle } \end{equation}
Turbulent Schmidt number over the entire flow depth. Dashed lines represent profiles from v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2025a ), where zero-flux scalar boundary conditions on the sphere surfaces were applied.
and its profiles are plotted in figure 20. In agreement with previous studies (e.g. Schwertfirm & Manhart Reference Schwertfirm and Manhart2007; Pirozzoli et al. Reference Pirozzoli, Bernardini and Orlandi2016),
${\textit{Sc}}_t \lesssim 1$
over large parts of the free-flow region. The different top domain boundary conditions for velocity and scalar leads to a rapid increase of
${\textit{Sc}}_t$
. Under adiabatic scalar boundary conditions on the surfaces of the porous medium, the turbulent Schmidt number
${\textit{Sc}}_t$
was found to decrease considerably in the near-interface region (e.g. Chandesris et al. Reference Chandesris, D’Hueppe, Mathieu, Jamet and Goyeau2013). The dashed lines in figure 20 represent profiles from our previous study (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
) and support this trend. In the present study with isothermal sphere surfaces, a slight, but critically less drastic, near-interface decrease of
${\textit{Sc}}_t$
is observed. As discussed in § 4.7, correlations between
$c^\prime$
and
$u^\prime$
could be seen (figure 19) and quantified (figure 18). As both fluctuations are likely to be connected to vertical velocity fluctuations
$w^\prime$
, the similar behaviour of the scalar transport
$\langle \overline {c^\prime w^\prime } \rangle$
and the turbulent shear stress
$\langle \overline {u^\prime w^\prime } \rangle$
can be explained. Due to the role of pressure drag, momentum absorption happens more effectively than scalar absorption. Particular for high
${\textit{Re}}_{{\kern-1pt}K}$
, this effect is reflected by different gradients
$\langle \partial \overline {c} \kern-0.5pt/\kern-1pt \partial z \rangle \lt \langle \partial \overline {u} \kern-0.5pt/\kern-1pt \partial z \rangle$
near the interface, which is likely to explain the slight near-interface decrease of
${\textit{Sc}}_t$
. The idea of different gradients is further supported by the spread of the respective absorption processes, which can be seen as a proxy for a characteristic length scale. The observation that
$\hat {\sigma }_z \gt \sigma _{z}$
(figure 11) implies that the scalar absorption happens more gradually than the momentum absorption, thus reducing the near-interface gradient. As the difference between
$\hat {\sigma }_z$
and
$\sigma _{z}$
appears to amplify with increasing
${\textit{Re}}_{{\kern-1pt}K}$
, the local decrease in
${\textit{Sc}}_t$
is also expected to become more prominent at higher
${\textit{Re}}_{{\kern-1pt}K}$
. As also remarked by Hantsis & Piomelli (Reference Hantsis and Piomelli2024),
${\textit{Sc}}_t$
has the drawback that it becomes ill-defined with increasing depth below the interface, as the denominator in (5.1) approaches zero.
Due to the no-slip boundary condition on the sphere surfaces, only diffusive scalar fluxes across the fluid–solid interface are possible (2.10), while advective scalar transport can additionally occur within the fluid-filled domain. According to the dimensionless form of the steady-state advection–diffusion equation (2.3), the generic Péclet number
${\textit{Pe}} = \textit{Re} \, \textit{Sc}$
determines the weighting between the advective and the diffusive term. Previously (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
), we argued that specifically the permeability Péclet number
${\textit{Pe}}_K = Re_K \, Sc$
is characteristic for the near-interface region, and determines the transition from the diffusive to the dispersive and further to the turbulent transport regime, as described by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2018). Accordingly, we assume that
${\textit{Pe}}_K$
has great influence on the region where the scalar is absorbed by the porous medium. While the present study considered different permeability Reynolds numbers
${\textit{Re}}_{{\kern-1pt}K} \in [0.4, 2.8]$
, we can only speculate about the impact of different Schmidt numbers. Whereas low values of
${\textit{Sc}} \ll 1$
are expected to enhance the near-surface scalar absorption, large values of
${\textit{Sc}} \gg 1$
at the same Reynolds number should allow the scalar flux to entrain deeper into the porous medium. From previous observations (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
), we expect that the dispersive transport would gain importance below the sediment bed surface, while affecting also deeper regions of the porous medium. It remains to be seen, however, if the categorisation of transport regimes suggested by Voermans et al. (Reference Voermans, Ghisalberti and Ivey2018) is finally recovered. As these ideas are speculative, an investigation of varying, and particularly of higher, Schmidt numbers would be an interesting extension of the present study.
6. Conclusion
In the present study, we investigated the absorption of a passive scalar on the surface of an isoscalar random sphere pack, which represents a porous bed underneath a turbulent free flow. An approach by single-domain pore-resolved DNS was chosen to simulate eight cases with different combinations of friction Reynolds numbers
${\textit{Re}}_\tau \in [150, 500]$
and permeability Reynolds numbers
${\textit{Re}}_{{\kern-1pt}K} \in [0.4, 2.8]$
. The varying ratios between flow depth and sphere diameter of
$h/D \in \{ 3, 5, 10 \}$
further allowed us to cover roughness Reynolds numbers of
$k_s^+ \in [20,200]$
, such that both transitionally and fully rough cases could be analysed. Based on the incompressible flow field, the advection–diffusion equation for a passive scalar with a unit Schmidt number
${\textit{Sc}} = 1$
was solved. The applied Dirichlet-type isoscalar boundary condition on the sphere surfaces contrasts a previous study (v. Wenczowski & Manhart Reference v. Wenczowski and Manhart2025a
), in which adiabatic sphere surfaces were defined and prevented scalar absorption.
The dimensionless scalar flux between the top boundary of the domain and the isoscalar sphere pack, expressed by the Sherwood number, increases primarily with
${\textit{Re}}_\tau$
. For cases with comparable
${\textit{Re}}_\tau$
, however, an evaluation of the Sherwood number showed that higher values of
${\textit{Re}}_{{\kern-1pt}K}$
enhance the scalar flux, hinting at more effective scalar absorption with increasingly permeable wall behaviour. Outer-layer similarities of the dimensionless velocity and scalar profiles allowed to determine the values of
$\Delta u^+$
and
$\Delta c^+$
, which represent roughness functions of the respective quantity. While
$\Delta c^+ \lt \Delta u^+$
is reported for rough walls, we observe that
$\Delta c^+ \approx \Delta u^+$
, which aligns with the enhanced scalar absorption due to permeability. In comparison to the absorption of momentum by the sphere pack (i.e. the action of drag), the scalar absorption peaks at a lower position and spreads over a wider vertical region. Both trends become more emphasised with increasing
${\textit{Re}}_{{\kern-1pt}K}$
, which can be explained by the role of pressure drag. From the mean position and the spread of the scalar absorption, we constructed a coordinate framework, in which similarities of the near-interface scalar field could be recovered.
Both no-slip velocity and isoscalar boundary conditions act as Dirichlet-type boundary conditions and prescribe uniform fixed values on the sphere surfaces, thus preparing the ground for several parallels between the scalar and velocity field: (i) absorption of scalar and momentum, respectively, are characteristic for the near-interface region, such that their analysis allows to construct coordinate frameworks, in which similarities can be recovered; (ii) within the sphere pack, spatial variations in both scalar and velocity field are moderate, as they can only develop on length scales associated with the pore space; (iii) steep bed-normal gradients prevail near the surface of the sphere pack, giving rise to gradient or shear production of temporal scalar or velocity fluctuations, respectively, while form-induced production plays a secondary role; (iv) near-interface temporal fluctuations of streamwise velocity and temporal fluctuations in the scalar field tend to be correlated and exhibit similar patterns at visual inspection. The similarity between scalar and velocity field is reflected by the turbulent Schmidt number
${\textit{Sc}}_t$
, which only shows a minor local decrease near the crests of the sphere pack.
Major differences prevail in comparison to scalar transport at
${\textit{Sc}} = 1$
, but with adiabatic conditions on the sphere surfaces. Under adiabatic scalar boundary conditions, regions of up- or down-welling fluid motion can induce strong spatial variations in the scalar field, such that dispersive scalar transport also affects deep layers of the sediment bed. By accelerating scalar diffusion, the isoscalar boundary condition drastically reduces spatial variations in the scalar field, thus restricting dispersive transport to a comparatively shallow layer below the sediment–water interface. We conclude that the scalar boundary condition on the sphere surfaces has strong impact on the sub-interface scalar transport and should not be neglected in models for hyporheic scalar transport. The influence of the Schmidt number remains as an interesting aspect to be addressed in future studies.
Acknowledgements
We gratefully acknowledge the correspondence with Joey Voermans and Marco Ghisalberti, as well as discussions with Yoshiyuki Sakai, Lukas Unglehrt and Rainer Helmig.
Funding
The authors gratefully acknowledge the financial support of the DFG under grant no. MA2062/15-1. Computing time was granted by the Leibniz Supercomputing Centre on SuperMUC-NG (project ‘pn68vo’).
Declaration of interests
The authors report no conflict of interest.
Author contributions
Simulations and postprocessing were performed by S. v. Wenczowski. Both authors contributed to reaching conclusions and writing the paper.
Appendix A. Momentum absorption
The absorption of momentum, which corresponds to the action of drag, is discussed in greater detail by v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2024), while the main aspects are summarised here for completeness. The motion of the fluid is described by the incompressible Navier–Stokes equation (see (2.1)), which is double-averaged according to the rules described (2.7). For the streamwise momentum, the superficially double-averaged equation can be rearranged to
\begin{align} \left ( f_{p,1}^s + f_{\nu ,1}^s \right ) = \frac {\partial }{\partial z} \left ( - \theta \!\left \langle \nu \frac {\partial \, \overline {u}}{\partial z} \right \rangle + \theta \!\left \langle \widetilde {\overline {u}} \, \widetilde {\overline {w}} \right \rangle + \theta \!\left \langle \overline {u^\prime w^\prime } \right \rangle \right ) - \theta \!\left \langle g_x \right \rangle \!. \end{align}
In (A1), the two terms on the left-hand side represent the superficially double-averaged pressure drag and viscous drag, respectively. Thus, an evaluation of the right-hand yields the complete drag. The individual drag components can be formulated as boundary terms or, alternatively, in terms of spatial averages and their spatial derivatives, as expressed by
\begin{align}&\qquad\qquad f_{p,1}^s = \frac {1}{\rho } \, {\textit{BT}}_x^{\,s}( -\overline {p} ) = \left \langle - \frac {1}{\rho } \frac {\partial p}{\partial x} \right \rangle ^s\!, \end{align}
\begin{align}& f_{\nu ,1}^s = {\textit{BT}}_{j}^{\,s} \left ( \nu \frac { \partial \overline {u} }{ \partial x_{\!j} } \right ) = \left \langle \frac {\partial }{\partial x_{\!j}} \left ( \nu \frac {\partial \overline {u}}{\partial x_{\!j}} \right ) \right \rangle ^s - \frac {\partial }{\partial z} \left \langle \nu \frac {\partial \overline {u}}{\partial z} \right \rangle ^s \!. \\[9pt] \nonumber \end{align}
As shown in the appendix of v. Wenczowski & Manhart (Reference v. Wenczowski and Manhart2025b
), the immersed boundary implementation allows a first-order accurate evaluation of the pressure drag via the relations in (A2). Figure 21 displays the contributions of the pressure drag to the total drag for case S-150 (
${\textit{Re}}_{{\kern-1pt}K} \approx 0.4$
) and for case M-500 (
${\textit{Re}}_{{\kern-1pt}K} \approx 2.8$
). The comparison of both cases highlights that the contribution of the pressure drag increases with
${\textit{Re}}_{{\kern-1pt}K}$
.
Contribution of the pressure drag to the total drag (a) for case S-150 (
${\textit{Re}}_{{\kern-1pt}K} \approx 0.4$
) and (b) for case M-500 (
${\textit{Re}}_{{\kern-1pt}K} \approx 2.8$
). The friction velocity
$u_\tau$
and the sphere diameter
$D$
are used for normalisation.
Open access
