1. Introduction
In turbulent channel flow, the wall plays an important role as a net source of vorticity injection into the domain. When considered from an observer’s perspective, the wall vorticity depends on the earlier-in-time flow state. Previous efforts have examined the contributions to the spanwise wall vorticity, or the streamwise surface stress, by stretching and tiling of earlier vorticity (Xiang, Eyink & Zaki Reference Xiang, Eyink and Zaki2025). In the present work, we are concerned with the utility of measurements of the wall stress to estimate the earlier flow state. Our objective is to determine the domain of dependence (DOD) of the measurement (Wang & Zaki Reference Wang and Zaki2025b ), which is defined as the earlier-in-time spatio-temporal state that can optimally alter the recorded surface stress. We develop an eddy-viscosity model to efficiently estimate the ensemble-averaged DOD in turbulent channel flow, compare it with the computationally expensive ensemble average, and demonstrate its universality and dynamics across Reynolds numbers.
The concept of a DOD originates from the theory of hyperbolic partial differential equations, where the solution at a given space–time point depends only on the initial condition within a characteristic cone with boundaries defined by the characteristic speeds of the system (Courant & Hilbert Reference Courant and Hilbert1962). For the Navier–Stokes equations, which are parabolic–elliptic and propagate information at infinite speed of sound, a sensitivity-based analogue is more appropriate: one seeks not which points affect the measurement, but the ones that have a strong causal influence. Identifying earlier flow events causally affecting a surface measurement is central to state estimation, model-based flow control and the physical interpretation of the dynamics. The adjoint framework is powerful for data assimilation and, closely related, the evaluation of the sensitivity of measurements to past flow events. In data assimilation, the adjoint field yields the gradient of the cost function with respect to the unknown parameters, and thus guides the minimisation of the mismatch between the simulation and measurements (Wang, Hasegawa & Zaki Reference Wang, Hasegawa and Zaki2019; Wang & Zaki Reference Wang and Zaki2025a ; Zaki Reference Zaki2025). In sensitivity analysis, the adjoint field is the flow state that has the greatest impact on an observable, and offers profound physical insight, particularly in linear, non-chaotic settings such as passive scalar transport in laminar flows (Pudykiewicz Reference Pudykiewicz1998; Neupauer & Wilson Reference Neupauer and Wilson2001; Keats, Yee & Lien Reference Keats, Yee and Lien2007; Wang & You Reference Wang and You2026). In turbulent flows, adjoint methods have demonstrated the inherent difficulty of reconstructing the full flow field from limited measurements. When the observations are the wall stresses, only near-wall structures and the largest outer-layer motions can be reliably recovered (Wang, Wang & Zaki Reference Wang, Wang and Zaki2022).
A key challenge in turbulent flows is the exponential amplification of the adjoint field in backward time, which is the dual to the classical butterfly effect (Deissler Reference Deissler1987; Lorenz Reference Lorenz2000). Infinitesimally close measurements can arise due to entirely different flow states (Zaki & Wang Reference Zaki and Wang2021). The exponential growth of each adjoint realisation in backward time can render both data assimilation and sensitivity analyses unstable or infeasible. The problem is especially acute when the target quantities are long-time averages (Garai & Murman Reference Garai and Murman2021; Yuan et al. Reference Yuan, Wang, Wang and Wang2023; Xu et al. Reference Xu, Zhao, Wu, Zhang, Müller, Huang, Rahmati and Wang2023). The difficulty has its roots in the fluctuation–response relations of statistical mechanics. Kubo’s theorem (Kubo Reference Kubo1957) states that the linear response of an equilibrium observable to a small perturbation equals an equilibrium time-correlation function of the unperturbed system. In the present context, the statistically stationary turbulent channel flow is the equilibrium ensemble, and Kubo’s theorem identifies the ensemble-averaged adjoint as a two-point space–time velocity correlation of the unperturbed turbulence. Two consequences follow. First, the mean adjoint must be smooth and well behaved even though individual realisations diverge exponentially, because the underlying equilibrium correlations of stationary turbulence are finite and regular. Second, any closure that accurately captures the dominant space–time correlations of the turbulence should yield a faithful approximation of the mean DOD without requiring direct ensemble averaging of chaotic trajectories.
To address these difficulties, stabilisation strategies have been developed. Least-squares shadowing (LSS) replaces the ill-conditioned initial-value problem with a well-posed variational problem, yielding accurate sensitivities of long-time averages in a simple, chaotic system (Wang Reference Wang2014; Wang, Hu & Blonigan Reference Wang, Hu and Blonigan2014; Chater et al. Reference Chater, Ni, Blonigan and Wang2017; Ni & Wang Reference Ni and Wang2017). The primal trajectory is replaced by a shadow trajectory that remains as close as possible to the original chaotic orbit while satisfying the linearised governing equations. The ill-posed backward-time initial-value problem is then converted into a coupled forward–backward boundary-value problem with a well-conditioned normal equation, eliminating the exponential divergence at the cost of a global time-window solve. Least-squares shadowing yields accurate sensitivities for low-dimensional chaotic systems such as the Lorenz equations and the Kuramoto–Sivashinsky equation (Wang et al. Reference Wang, Hu and Blonigan2014; Chater et al. Reference Chater, Ni, Blonigan and Wang2017), but its computational cost scales with both the state dimension and the time horizon, making direct application to turbulent flows prohibitive. Non-intrusive least-squares shadowing (Ni & Wang Reference Ni and Wang2017) reduces this cost by restricting the shadow correction to the unstable subspace of the flow, requiring only a number of additional forward simulations equal to the number of positive Lyapunov exponents. However, the dimension of the unstable subspace grows with Reynolds number in turbulent shear flows, limiting the practical usage (Chandramoorthy et al. Reference Chandramoorthy, Fernandez, Talnikar and Wang2019). Adjoint-energy-based stabilisation introduces artificial damping motivated by the adjoint energy budget, and offers a cheaper alternative for turbulent flows (Blonigan et al. Reference Blonigan, Chen, Wang and Larsson2012; Garai & Murman Reference Garai and Murman2021).
In the context of the present work, the mean DOD requires an ensemble average of the adjoint fields (Wang et al. Reference Wang, Wang and Zaki2022), evaluated from different realisations of the forward trajectory. However, due to a power-law tail in the histogram of adjoint gradients, ensemble averaging converges very slowly (Eyink, Haine & Lea Reference Eyink, Haine and Lea2004). A potential remedy to this slow convergence is sought, and is inspired by the literature on the forward-time linearised dynamics in turbulent channel flow. Rather than study the forward evolution of perturbations to the full turbulent state, which are exponentially unstable, Del Alamo & Jimenez (Reference Del Alamo and Jimenez2006) and Cossu, Pujals & Depardon (Reference Cossu, Pujals and Depardon2009) considered the mean turbulent profile and adopted an eddy-viscosity model for the analysis. They demonstrated that, despite the damping of all eigenmodes, perturbations can exhibit significant transient energy growth that closely mirrors the organisation and scaling of coherent structures observed in wall turbulence. Vadarevu et al. (Reference Vadarevu, Symon, Illingworth and Marusic2019) showed that an impulsive body force in the linearised Navier–Stokes system, augmented with an eddy-viscosity model, produces vortex–streak structures that exhibit geometric self-similarity and decay predictably with time. These analyses predict an inner and outer peak in the optimal energy amplification, and a self-similar evolution of streamwise streaks. In addition to the growth of streaks by lift-up (Landahl Reference Landahl1980; Liu et al. Reference Liu, Semin, Godoy-Diana and Wesfreid2024b ), Orr amplification is another important non-modal mechanism where perturbations that are initially tilted against the mean shear grow as they are rotated toward the direction of maximum strain. Encinar & Jiménez (Reference Encinar and Jiménez2020) demonstrated that conditionally averaged burst structures in fully nonlinear direct numerical simulations (DNS) are consistent with the Orr reorientation of upstream-tilted perturbations. Their results establish a direct connection between observed turbulent structures and the non-normal amplification of linearised dynamics. While the above studies focus on the forward-time evolution and statistical interpretation of coherent structures, the present work addresses the complementary problem of backward-time sensitivity: namely, which earlier perturbations most strongly influence a given wall-stress measurement.
The rest of the paper is structured as follows. In § 2, we present the mathematical formulation for the eddy-viscosity closure of the ensemble-averaged linear forward and adjoint equations, and provide a physical interpretation of the mean adjoint fields. In § 3, we contrast the forward and adjoint fields, and demonstrate their respective self-similar behaviours. We compare the adjoint model with ensemble averaging of realisations, and examine the forward-time dynamics of the computed adjoint DOD to aid the interpretation. Conclusions are provided in § 4.
2. Mathematical formulation
We start from the linearised incompressible Navier–Stokes equations for turbulent channel flow, and perform Reynolds averaging to evaluate the ensemble statistics. This leads to the Reynolds-averaged linearised Navier–Stokes equations (RA-LNS). With an appropriate closure model, we obtain the forward formulation and then derive the associated Reynolds-averaged adjoint Navier–Stokes (RA-ANS) system for backward-time sensitivity analysis. The reasoning and limitations of eddy-viscosity closure are discussed. The physical interpretations of these equations are then presented: the RA-LNS describes the forward-in-time domain of influence (DOI) of a disturbance, while the RA-ANS quantifies the back-in-time DOD of an observation. All the numerical solutions of the forward and adjoint Reynolds-averaged equations reported herein are performed in spectral space, using standard Fourier–Fourier–Chebyshev discretisation of the governing equations.
2.1. The Reynolds-averaged linearised Navier–Stokes equations
Suppose the incompressible velocity field
$\boldsymbol{U}(\boldsymbol{x},t)$
is a realisation of statistically stationary turbulent channel flow. Variables are non-dimensionalised by the half-channel height
$h^\star$
, the ensemble-averaged bulk velocity
$\left \langle U_b^\star \right \rangle$
and the fluid density
$\rho ^\star$
, where superscript
$\star$
marks dimensional quantities and
$\langle \bullet \rangle$
is the ensemble average defined below. The non-dimensional Navier–Stokes equations are
\begin{eqnarray} {\partial _t \boldsymbol{U}} + \left (\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla } \right ) \boldsymbol{U} &=& - \boldsymbol{\nabla } P + \nu {\nabla} ^2 \boldsymbol{U}, \qquad \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{U} = \boldsymbol{0}, \end{eqnarray}
where the non-dimensional kinematic viscosity
$\nu$
and bulk Reynolds number are related by
$Re = 1/\nu := \left \langle U^{\star }_b \right \rangle h^{\star } / \nu ^{\star }$
. The non-dimensional kinematic pressure is
$P = P^\star/\rho^\star U_b^{\star 2}$
. Superscript
$+$
denotes normalisation by viscous scales, specifically, the ensemble-averaged friction velocity
$\left \langle U_{\tau }^\star \right \rangle := \sqrt { \nu ^\star ({\partial \left \langle U^\star \right \rangle }/{\partial y^\star } )_{y^\star =0} }$
and length scale
$\nu ^{\star }/\left \langle U_{\tau }^\star \right \rangle$
. The friction Reynolds number is defined as
$Re_\tau :=\left \langle U_\tau ^\star \right \rangle h^\star / \nu ^\star$
.
We introduce a small perturbation field
$(\boldsymbol{u},p)$
superimposed on top of
$(\boldsymbol{U},P)$
. This perturbation is thus governed by the linearised Navier–Stokes equations
\begin{equation} {\partial _t \boldsymbol{u}} + \left (\boldsymbol{U} \boldsymbol{\cdot }\boldsymbol{\nabla } \right ) \boldsymbol{u} + \left (\boldsymbol{u} \boldsymbol{\cdot }\boldsymbol{\nabla } \right ) \boldsymbol{U} = - \boldsymbol{\nabla } p + \nu{\nabla} ^2 \boldsymbol{u},\qquad \boldsymbol{\nabla } \boldsymbol{\cdot }\boldsymbol{u} = \boldsymbol{0},\qquad \boldsymbol{u}_0 = \boldsymbol{\phi }, \end{equation}
where
$\boldsymbol{\phi }$
is the initial velocity perturbation. The ensemble average of the above equations over different realisations of the turbulent field
$\boldsymbol{U}$
yields
\begin{align} &{\partial _t \left \langle \boldsymbol{u} \right \rangle } + \left \langle \boldsymbol{U} \right \rangle \boldsymbol{\cdot }\boldsymbol{\nabla } \left \langle \boldsymbol{u} \right \rangle + \left (\left \langle \boldsymbol{u} \right \rangle \boldsymbol{\cdot }\boldsymbol{\nabla } \right ) \left \langle \boldsymbol{U} \right \rangle = -\boldsymbol{\nabla }\boldsymbol{\cdot }\left (\left \langle \boldsymbol{U}^\prime \boldsymbol{u}^{\prime } +\boldsymbol{u}^{\prime }\boldsymbol{U}^\prime \right \rangle \right ) - \boldsymbol{\nabla } \left \langle p \right \rangle + \nu {\nabla} ^2 \left \langle \boldsymbol{u} \right \rangle, \end{align}
\begin{align} &\quad\quad\kern-5pt\boldsymbol{\nabla } \boldsymbol{\cdot }\left \langle \boldsymbol{u} \right \rangle = \boldsymbol{0},\qquad \left \langle \boldsymbol{u} \right \rangle _0 = \boldsymbol{\phi }. \end{align}
The ensemble average
$\langle \bullet \rangle$
is computed using realisations of statistically stationary turbulent base flows, and exploiting ergodicity. Specifically, a long trajectory of the forward flow is sampled evenly in time, with circular translation in the homogeneous
$x$
and
$z$
directions, to construct an ensemble of base flows
$\boldsymbol{U}(\boldsymbol{x},t)$
. Each member of this ensemble is then adopted as the base flow to advance the perturbation fields
$\boldsymbol{u}$
. In (2.3), the terms
$\boldsymbol{U}'\boldsymbol{u}'+\boldsymbol{u}'\boldsymbol{U}'$
are outer (dyadic) products, with components
$(U_i'u_j'+u_i'U_j')$
. The primed variable
$\boldsymbol{U}^\prime \equiv \boldsymbol{U} - \left \langle \boldsymbol{U} \right \rangle$
is the deviation away from the mean turbulent velocity field, and
$\boldsymbol{u}^\prime \equiv \boldsymbol{u} - \left \langle \boldsymbol{u} \right \rangle$
is the deviation from the mean perturbation field. Both deviations satisfy the divergence-free condition. We note from (2.2) and (2.3) the distinct dynamics of the evolution of an individual perturbation
$\boldsymbol{u}$
and of the mean
$\left \langle \boldsymbol{u} \right \rangle$
evaluated over different realisations of
$\boldsymbol{U}$
. As shown in previous studies (Nikitin Reference Nikitin2018; Wang et al. Reference Wang, Wang and Zaki2022), the instantaneous perturbation field grows exponentially in forward time at the Lyapunov rate. If we take a perturbed initial condition
$\boldsymbol{U}_0^{\delta } = \boldsymbol{U}_0 + \boldsymbol{u}_0$
, the evolution of this new field would gradually converge to the same statistically stationary channel-flow turbulence with the same mean, i.e.
$\left \langle \boldsymbol{U}^{\delta } \right \rangle =\left \langle \boldsymbol{U} \right \rangle$
as time
$t \rightarrow \infty$
, or equivalently
$\left \langle \boldsymbol{U}^\prime \right \rangle = 0$
. As such, the ensemble-averaged difference between the perturbed and the original fields must decay at long forward times. As for the linear perturbations, each realisation
$\boldsymbol{u}$
amplifies exponentially while the ensemble average
$\left \langle \boldsymbol{u} \right \rangle (t)$
decays in time (Malkus Reference Malkus1956; Markeviciute & Kerswell Reference Markeviciute and Kerswell2023; Zhu, Chen & Fu Reference Zhu, Chen and Fu2024). This distinction between fluctuation growth and mean decay underlies the difficulty of computing the ensemble average in order to evaluate the mean of the linearised fields, which requires an exponentially growing number of samples as the time of interest is increased.
We proceed to introduce a linear eddy-viscosity model (LEVM) to close
$\left \langle \boldsymbol{U}^\prime \boldsymbol{u}^{\prime } +\boldsymbol{u}^{\prime }\boldsymbol{U}^\prime \right \rangle$
. In the original LEVM, the turbulent stress tensor
$\left \langle \boldsymbol{U}^\prime \boldsymbol{U}^\prime \right \rangle$
is modelled as
\begin{equation} -\left \langle U_i^\prime U_j^\prime \right \rangle = 2\nu _t \left \langle S \right \rangle _{ij} - \frac 23 \left \langle K \right \rangle \delta _{ij}, \end{equation}
where
$\nu _t$
is the effective turbulent viscosity,
$\left \langle S \right \rangle _{ij} = ( 1/2) (\partial _j \left \langle U \right \rangle _i + \partial _i \left \langle U \right \rangle _j )$
is the mean strain-rate tensor for the base state
$\boldsymbol{U}$
and
$\left \langle K \right \rangle =( 1/2) \left \langle U_i^{\prime }U_i^{\prime } \right \rangle$
is the fluctuation kinetic energy in the base state, and
$\delta_{\textit{ij}}$
is the Kronecker delta. Here,
$\epsilon \ll 1$
is an order parameter characterising the magnitude of the velocity perturbation
$\boldsymbol{u}$
relative to the base-flow scale, i.e.
$|\boldsymbol{u}|=O(\epsilon |\boldsymbol{U}|)$
; the same
$\boldsymbol{u}$
appearing in ((2.2)–(2.3)) is used throughout. Introducing this perturbation into the eddy-viscosity relation yields the same statistical pattern at leading order in
$\epsilon$
, Adding a small perturbation of magnitude
$o(\epsilon )$
to this flow field yields the same statistical pattern
\begin{equation} -\left \langle (U_i^\prime + \epsilon u_i) (U_j^\prime + \epsilon u_j) \right \rangle = 2\nu _t (\left \langle S \right \rangle _{ij} + \epsilon \left \langle s \right \rangle _{ij}) - \frac 23 (\left \langle K \right \rangle + \epsilon \left \langle U_l^\prime u_l \right \rangle + \epsilon ^2 \left \langle k \right \rangle )\delta _{ij}. \end{equation}
Here,
$\left \langle s \right \rangle _{ij} = (1/2) (\partial _j \left \langle u \right \rangle _i + \partial _i \left \langle u \right \rangle _j )$
is the mean strain-rate tensor for the perturbation field and
$\left \langle k \right \rangle =(1/2) \left \langle u_i^{\prime }u_i^{\prime } \right \rangle$
is the fluctuation kinetic energy for the perturbation field
$\boldsymbol{u}$
. Subtracting (2.4) and (2.5), and dividing by
$\epsilon$
, yields
\begin{equation} -\left \langle U_i^\prime u_j +U_j^\prime u_i \right \rangle = 2\nu _t \left \langle s \right \rangle _{ij}- \frac 23 \left \langle U_l^\prime u_l \right \rangle \delta _{ij}+\mathcal{O}(\epsilon ). \end{equation}
For clarity, we note that the repeated
$l$
index implies summation (Einstein convention). Ignoring the higher-order terms and applying
$\partial _x \left \langle U \right \rangle = \partial _z \left \langle U \right \rangle = \left \langle V \right \rangle = \left \langle W \right \rangle = 0$
in a turbulent channel flow leads to equations for the modelled mean base flow
$\overline {U}$
and mean perturbations
$\overline {\boldsymbol{u}}$
and
$\overline {p}$
, as approximations of
$\left \langle U \right \rangle$
,
$\left \langle \boldsymbol{u} \right \rangle$
and
$\left \langle p \right \rangle$
\begin{equation} \partial _t \overline {\boldsymbol{u}} + \overline {U} \partial _x \overline {\boldsymbol{u}} + \overline {v} \;d_y \overline {U} \, \boldsymbol{e}_x = - \boldsymbol{\nabla } \overline {p} + \nu _t {\nabla} ^2 \overline {\boldsymbol{u}} + \left (d_y\nu _t\right ) \left (\partial _y \overline {\boldsymbol{u}} + \boldsymbol{\nabla }\overline {v}\right ), \quad \boldsymbol{\nabla } \boldsymbol{\cdot }\overline {\boldsymbol{u}} = 0, \quad \overline {\boldsymbol{u}}_0 = \boldsymbol{\phi }. \end{equation}
In the above expression, the pressure
$\overline {p}$
includes the
$\left \langle U_l^\prime u_l \right \rangle$
term. The velocity vector is
$\bar{\boldsymbol{u}}=(\bar{u},\bar{v},\bar{w})$
, where
$\bar{v}$
is the wall-normal component, and
$\boldsymbol{e}_x$
denotes the unit vector in the streamwise direction. This set of equations can be termed the RA-LNS, and they have the same form as the linearised Reynolds-averaged Navier–Stokes equations. In other words, the ensemble average of linear perturbations to different realisations of the base flow is the same as the linear perturbation to the ensemble-averaged base flow. The closure above implicitly assumes that the correlation term
$\left \langle \boldsymbol{U}'\boldsymbol{u}'+\boldsymbol{u}'\boldsymbol{U}' \right \rangle$
depends instantaneously on the mean perturbation strain
$\partial _y \overline {\boldsymbol{u}} + \boldsymbol{\nabla }\overline {v}$
through an eddy viscosity (Hamba Reference Hamba1995, Reference Hamba2004; Liu et al. Reference Liu, Schäfer, Bryngelson, Zaki and Mani2024a
). In principle, the interaction between the turbulent fluctuations
$\boldsymbol{U}'$
and the perturbation field
$\boldsymbol{u}$
may involve a non-local dependence in time, so that the effective turbulent stress would contain a memory kernel that depends on the past history of the perturbation. In the present work, we neglect such memory effects and adopt a steady eddy-viscosity closure that depends only on the instantaneous mean strain.
The herein adopted eddy-viscosity model is an approximation to the true Reynolds stress dynamics, and linearisation about this model introduces an additional simplification. Nonetheless, this approach has been successful in capturing the linear dynamics, including transient growth, and in explaining the organisation and scaling of coherent structures in turbulent flows (Del Alamo & Jimenez Reference Del Alamo and Jimenez2006; Cossu et al. Reference Cossu, Pujals and Depardon2009; Vadarevu et al. Reference Vadarevu, Symon, Illingworth and Marusic2019). For the present purposes, we model the eddy viscosity using the analytic expression by Cess (Reference Cess1958)
\begin{align} \frac {\nu _t(y)}{\nu }= \frac {1}{2}\left \{1+\frac {\kappa ^2 {Re}_\tau ^2}{9}\big (2y-y^2\big )^2\big (3-4y+2y^2\big )^2\right.\\\nonumber \quad\left.\left [1-\exp \left [(|1-y|-1) {Re}_\tau / A\right ]\right ] \vphantom{\frac {\kappa ^2 {Re}_\tau ^2}{9}}\right\}^{1 / 2}+\frac {1}{2}. \end{align}
The von Kármán constant is
$\kappa =0.426$
and
$A=25.4$
. The computation of the turbulent mean velocity profile
$\overline {U}(y)$
involves integration of
$(y-1)\overline {U}_{\tau } Re_{\tau }/(\nu _t/\nu )$
from
$0$
to
$y$
.
2.2. The Reynolds-averaged adjoint equations and the mean DOD
The solution of the averaged linearised forward (2.7) from the previous section can be written as the linear map
\begin{equation} \left \langle \boldsymbol{u} \right \rangle _{t_m} = \mathcal{A}\left \langle \boldsymbol{u} \right \rangle _0, \end{equation}
where
$\mathcal{A}$
is a linear operator acting on the space of divergence-free vector fields with sufficient smoothness, representing the evolution of the perturbation from the initial condition
$\left \langle \boldsymbol{u} \right \rangle _0$
to the measurement time
$t_m$
. The adjoint evolution
$\mathcal{A}^{\dagger }$
is then derived from the forward–adjoint duality relation
\begin{equation} \mathcal{I}\!\left ( \mathcal{A}\left \langle \boldsymbol{u} \right \rangle _0, \left \langle \boldsymbol{u} \right \rangle ^{\dagger }_{t_m} \right ) = \mathcal{I}\!\left ( \left \langle \boldsymbol{u} \right \rangle _0, \mathcal{A}^{\dagger }\left \langle \boldsymbol{u} \right \rangle ^{\dagger }_{t_m} \right ), \end{equation}
where the inner product is defined as
$\mathcal{I}\! (\boldsymbol{f}, \boldsymbol{g} ) := \int _V \boldsymbol{f}^H\boldsymbol{g} {\rm d}V$
and
$\left \langle \boldsymbol{u} \right \rangle ^{\dagger }$
is the adjoint of the ensemble-averaged velocity perturbation. The RA-ANS equations for the modelled mean perturbation fields
$\overline {\boldsymbol{u}}$
, are derived from this duality relation using integration by parts, and are given by
\begin{align} &\partial _{\zeta } \overline {\boldsymbol{u}}^{\dagger } - \overline {U} \partial _x \overline {\boldsymbol{u}}^{\dagger } + d_y\overline {U} \overline {u}^{\dagger } \boldsymbol{e}_y = \boldsymbol{\nabla } \overline {p}^{\dagger } + \nu _t {\nabla} ^2 \overline {\boldsymbol{u}}^{\dagger } + \left (d_y\nu _t\right ) \partial _y \overline {\boldsymbol{u}}^{\dagger } - \left (d_y^2\nu _t\right ) \overline {v}^{\dagger } \boldsymbol{e}_y, \end{align}
\begin{align} & \boldsymbol{\nabla } \boldsymbol{\cdot }\overline {\boldsymbol{u}}^{\dagger } = \boldsymbol{0}, \qquad \overline {\boldsymbol{u}}^{\dagger }_{t_m} = \boldsymbol{\phi }^{\dagger }. \end{align}
The first equation is thus the eddy-viscosity modelled version of the ensemble-averaged adjoint momentum equation
\begin{equation} \partial _{\zeta } \left \langle \boldsymbol{u} \right \rangle ^{\dagger } - \left \langle \boldsymbol{U} \right \rangle \boldsymbol{\cdot }\boldsymbol{\nabla }\left \langle \boldsymbol{u} \right \rangle ^{\dagger } + \boldsymbol{\nabla }\left \langle \boldsymbol{U} \right \rangle \boldsymbol{\cdot }\left \langle \boldsymbol{u} \right \rangle ^{\dagger } = \boldsymbol{\nabla } \left \langle p \right \rangle ^{\dagger } + \nu {\nabla} ^2 \left \langle \boldsymbol{u} \right \rangle ^{\dagger } + \left \langle \boldsymbol{U}^{\prime }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^{\dagger } - \boldsymbol{\nabla }\boldsymbol{U}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\dagger } \right \rangle . \end{equation}
The adjoint equations are solved in backward time,
$\zeta = t_m-t$
. The appearance of the transposed advection operator
$\boldsymbol{\nabla }\left \langle \boldsymbol{U} \right \rangle \boldsymbol{\cdot }\left \langle \boldsymbol{u} \right \rangle ^{\dagger }$
fundamentally distinguishes the adjoint dynamics from the forward system (Wang et al. Reference Wang, Wang and Zaki2022). In addition, the correlation term
$\left \langle \boldsymbol{U}^{\prime }\boldsymbol{\cdot }\boldsymbol{\nabla }\boldsymbol{u}^{\dagger } - \boldsymbol{\nabla }\boldsymbol{U}^{\prime }\boldsymbol{\cdot }\boldsymbol{u}^{\dagger } \right \rangle$
represents the interaction between forward fluctuations and the adjoint field. Such forward–adjoint correlations are rarely considered in turbulence modelling, and their statistical closure remains largely unexplored. A comparison of (2.7) and (2.11) reveals that the eddy-viscosity closure is not self-adjoint in the presence of mean shear. The forward eddy-viscosity term in (2.7) is
$\nu _t{\nabla} ^2\overline {\boldsymbol{u}} + (d_y\nu _t )(\partial _y\overline {\boldsymbol{u}}+\boldsymbol{\nabla }\overline {v})$
, whereas the corresponding term in the adjoint (2.11) is
$\nu _t{\nabla} ^2\overline {\boldsymbol{u}}^{\dagger } + (d_y\nu _t )\,\partial _y\overline {\boldsymbol{u}}^{\dagger } - (d_y^2\nu _t )\,\overline {v}^{\dagger }\boldsymbol{e}_y$
. The extra term
$- (d_y^2\nu _t )\,\overline {v}^{\dagger }\boldsymbol{e}_y$
arises from integration by parts of the wall-normal gradient contribution and reflects the wall-normal inhomogeneity of
$\nu _t(y)$
. In the limiting case of a spatially uniform eddy viscosity (
$d_y\nu _t = 0$
), both operators reduce to
$\nu _t{\nabla} ^2$
, which is self-adjoint; it is the wall-normal variation of
$\nu _t$
driven by mean shear that breaks the self-adjoint structure of the closure.
To demonstrate the physical meaning of the adjoint field, we prescribe its terminal condition
$\boldsymbol{u}^{\dagger }_{t_m} = \boldsymbol{\phi }^{\dagger }$
as the measurement kernel, i.e.
$\mathcal{I}\! (\left \langle \boldsymbol{u} \right \rangle _{t_m}, \boldsymbol{\phi }^{\dagger } ) = \left \langle m \right \rangle$
, and substitute into the duality relation. We obtain
\begin{equation} \mathcal{I}\!\left (\left \langle \boldsymbol{u} \right \rangle _{t_m}, \right ){\boldsymbol{\phi }^{\dagger } }= \mathcal{I}\!\left (\mathcal{A}\left \langle \boldsymbol{u} \right \rangle _0, \boldsymbol{\phi }^{\dagger } \right ) = \mathcal{I}\!\left (\left \langle \boldsymbol{u} \right \rangle _0, \mathcal{A}^{\dagger }\boldsymbol{\phi }^{\dagger } \right ) = \mathcal{I}\!\left (\left \langle \boldsymbol{u} \right \rangle _0, \left \langle \boldsymbol{u} \right \rangle _0^{\dagger } \right ). \end{equation}
As such, the mean variation in the measurement,
$\left \langle m \right \rangle$
, is proportional to the mean initial perturbation
$\left \langle \boldsymbol{u} \right \rangle _0$
, and the proportionality is given by the mean adjoint field
$\left \langle \boldsymbol{u} \right \rangle _0^{\dagger }$
. The mean adjoint can therefore be interpreted as the mean DOD for the measurement kernel
$\boldsymbol{\phi }^{\dagger }$
. This idea is demonstrated schematically in figure 1, where the DOI and DOD using one realisation of channel-flow turbulence are plotted, along with their ensemble averages, which are shown in transparent iso-surfaces.
Examples of the instantaneous and mean DOI from the linearised Navier–Stokes equations and DOD from the adjoint equations. Opaque isosurfaces: DOI and DOD computed from one realisation of channel-flow turbulence. Transparent isosurfaces: mean DOI and DOD.
3. Results
3.1. The mean domain of dependence
We solve both the Reynolds-averaged linearised and adjoint Navier–Stokes (ANS) equations, RA-LNS (2.7) and RA-ANS (2.11), respectively, initialised using wall-stress data. The forward equations provide the mean DOI of a wall-stress perturbation, and the adjoint equations provide the mean DOD of a wall-stress measurement kernel in backward time. The base state for both the RA-LNS and RA-ANS equations is the mean base-flow profile
$\overline {U}(y)$
, computed using the Cess eddy-viscosity model (2.8). This profile satisfies the steady Reynolds-averaged Navier–Stokes equations for fully developed and statistically stationary channel flow, and not the instantaneous Navier–Stokes equations (2.1). The initial wall-stress kernel at the wall is the derivative of a narrow Gaussian that represents
$\partial u/\partial y|_{y=0}$
. Specifically, we use
$ \phi (\boldsymbol{x}) = ({y^+}/{\sigma ^5 ( 2 \pi )^{3/2}} ) \exp (-{|\boldsymbol{x}|^2}/{2 \sigma ^2} )$
, where
$\sigma$
is defined as the width of the Gaussian in viscous units. A projection operator that enforces the solenoidal and boundary conditions is applied to the initial condition. The forward solver is validated against the results from Pujals et al. (Reference Pujals, García-Villalba, Cossu and Depardon2009), including optimal transient growth curves. The adjoint solver is validated using the forward–adjoint duality relation described in § 2. The friction Reynolds numbers that were examined are
$Re_{\tau } = \{180,590,10^3,10^4\}$
. For a given friction Reynolds number, the linearised and adjoint fields are evaluated for streamwise wavenumbers
$k_x^+ \in [0,0.4775]$
sampled at intervals of
$\varDelta k^+_x=2.5\times 10^{-3}$
, and spanwise wavenumbers
$k_z^+ \in [-0.5,0.5]$
sampled at intervals of
$\varDelta k_z^+ = 4\times 10^{-3}$
.
Comparison of the mean DOI and DOD at
$Re_\tau = 10^3$
. (a) The DOI (left column) and the DOD (right column) starting from the measurement kernel of
$ {\partial u}/{\partial y}$
at the wall (
$y=0$
). Iso-surfaces of the streamwise velocity components of these structures are shown at selected forward- and backward-time instances. The scale in the upper-left corner applies to all panels. (b) Contours of the energy in different wavenumbers
$(k_x^+, k_z^+)$
, and vertically integrated across the channel, are shown for the DOI (left column) and DOD (right column) at the corresponding instances.
Iso-surfaces of the normalised mean DOI and DOD are shown in figure 2
$(a)$
, and their energy spectra are reported in panel
$(b)$
. The results are for
$Re_{\tau } = 10^3$
, and those for the other Reynolds numbers are visually indistinguishable when plotted in viscous units. The DOI and DOD expand in forward and reverse time, respectively, while largely preserving their shapes (figure 2
a). The qualitative difference in their form in physical space reinforces that the forward influence of a perturbation of wall stress and the sensitivity of a wall-stress measurement are not mirror images of each other. In wavenumber space (panel
$(b)$
), the energy distribution of the DOI is concentrated in a single, simply connected region, characterised by a dominant wavenumber. In contrast, the DOD energy is distributed across two distinct regions: one region is associated with upstream-leaning
$u$
-velocity sheets, which correspond to low
$k_z$
in the spectra. This region is dominant at early backward times and persists at longer backward times. The second region is associated with streamwise streaky structures. These are low
$k_x$
in the spectra, and become appreciable at longer backward times. The dynamics of these two regions will be examined after a comparison of the mean DOD and the ensemble average from adjoint realisations.
Comparison of the isosurfaces (a)
$\overline {u}^{\dagger }/\overline {u}^{\dagger }_{max} = \pm 0.01$
and (b)
$\langle u^\dagger \rangle /\langle u^\dagger \rangle _{max} = \pm 0.01$
, computed from the RA-ANS and the ensemble adjoint approach (5000 samples), respectively, for backward times
$\zeta ^+ = \{4, 8, 12, 16, 20, 50\}$
for
$Re_{\tau }=180$
. The background grey colour shows a sample of the underlying turbulent streamwise wall-shear stress.
Figure 3. Long description
Panel A: Three-dimensional heat maps showing isosurfaces computed from the RA-ANS eddy-viscosity model. The color scale ranges from -0.01 to 0.01, with blue representing lower values and red representing higher values. The isosurfaces are displayed at different backward times, with the background showing the turbulent streamwise wall-shear stress in grey. Panel B: Three-dimensional heat maps showing isosurfaces computed from the ensemble adjoint approach with 5000 samples. Similar to Panel A, the color scale ranges from -0.01 to 0.01, with blue representing lower values and red representing higher values. The isosurfaces are displayed at different backward times, with the background showing the turbulent streamwise wall-shear stress in grey.
(a.i–iv) Instances of the DOD, shown by isosurfaces of
${u}^{\dagger }/{u}^{\dagger }_{max} = \pm 0.01$
, for four different realisations of the turbulent base flow, at the backward times
$\zeta ^+ = \{4, 8, 12, 16, 20, 50\}$
for
$Re_{\tau }=180$
. (b) Averaged adjoint structure using ten ensemble members.
To verify the eddy-viscosity adjoint model, we compare its prediction of the mean adjoint field with an ensemble average of 5000 realisations from ANS simulations. For each ensemble member, the adjoint computation adopts an independent realisation of the turbulent base flow
$\boldsymbol{U} (\boldsymbol{x},t )$
, sampled from a long statistically stationary DNS of channel flow. The base states were sampled at fifty instants separated by
$\varDelta t^+=20$
, and the fields were circularly translated ten times in each of the horizontal
$x$
and
$z$
directions. Figure 3 shows the isosurfaces of the streamwise adjoint velocity from both approaches. The structures exhibit strong visual similarity, particularly for short backward times (
$\zeta ^+ \le 20$
). At longer times, such as
$\zeta ^+ = 50$
, the ensemble-averaged adjoint field features high-frequency oscillations which are indicative of a lack of statistical convergence.
In figure 4(
$a$
), we report four realisations of the adjoint from the ANS simulations, which should be viewed in contrast to the ensemble average in figure 3. The individual realisations exhibit appreciable variability among themselves, with complex small-scale structures due to the chaotic backward-time evolution. These small scales are, however, absent from the ensemble average in figure 3, which captures the coherent mean structure. Figure 4(
$b$
) is the ensemble average of ten realisations only, which should again be viewed with reference to the converged average using 5000 samples in figure 3. The comparison underscores that the statistical convergence of the ensemble average becomes progressively more challenging at longer backward time.
In order to further quantify the alignment, or similarity, between adjoint fields, we define the spatial correlation coefficient
\begin{equation} \mathcal{R}(\boldsymbol{a},\boldsymbol{b}) =\frac {\int _V \boldsymbol{a} \boldsymbol{\cdot }\boldsymbol{b}\,{\rm d}V} {\left (\int _V|\boldsymbol{a}|^2\,{\rm d}V\right )^{1/2} \left (\int _V|\boldsymbol{b}|^2\,{\rm d}V\right )^{1/2}}. \end{equation}
In figure 5, we report the correlation between the RA-ANS
$\bar {\boldsymbol{u}}^{\dagger }$
and ten adjoint realisations
$\boldsymbol{u}^{\dagger }$
(grey curves). These correlations are weak, often of the order of
$0.1$
and at most approximately
$0.5$
, which is due to the chaotic backward-time evolution of individual adjoint realisations. The figure also reports the correlation between the RA-ANS
$\bar {\boldsymbol{u}}^{\dagger }$
and the ensemble average
$\langle \boldsymbol{u}^{\dagger } \rangle$
from
$5000$
realisations (black line). This correlation remains strongly positive, reaching approximately
$0.7$
over the time interval where the ensemble statistics are converged. This result demonstrates that while individual adjoint trajectories are highly variable, the RA-ANS formulation successfully captures the dominant direction of the ensemble-averaged sensitivity.
Spatial correlation coefficient between the RA-ANS adjoint field
$\bar {\boldsymbol{u}}^{\dagger }$
and the ensemble-averaged adjoint
$\langle \boldsymbol{u}^{\dagger } \rangle$
obtained from
$5000$
realisations (black line), and ten representative individual adjoint realisations
$\boldsymbol{u}^{\dagger }_{(n)}$
without averaging (grey lines). The strong positive correlation (approximately
$0.7$
) with the ensemble mean, contrasted with the weak correlation with individual realisations, demonstrates that RA-ANS captures the ensemble-averaged sensitivity direction while filtering the chaotic fluctuations of single adjoint trajectories.
The alignment between the RA-ANS
$\bar {\boldsymbol{u}}^{\dagger }$
and the ensemble average
$\langle {\boldsymbol{u}}^{\dagger }\rangle$
has implications for adjoint-variational data assimilation, where the gradient of the cost functional is furnished by the adjoint field. In that context, the adjoint field establishes the direction of steepest ascent (or descent) of the cost functional. Although individual adjoint realisations exhibit chaotic divergence and are poorly correlated with the RA-ANS field, the strong positive correlation between
$\bar {\boldsymbol{u}}^{\dagger }$
and the ensemble-averaged adjoint
$\langle \boldsymbol{u}^{\dagger } \rangle$
indicates that the RA-ANS solution captures the dominant gradient direction associated with the mean sensitivity. As a result, updates proportional to
$\bar {\boldsymbol{u}}^{\dagger }$
are expected to produce a positive ensemble-averaged first-order variation of the measurement functional. In this sense, the RA-ANS formulation can be viewed as a regularised approximation of the ensemble sensitivity direction, analogous to a regularised gradient (Protas, Bewley & Hagen Reference Protas, Bewley and Hagen2004), or a preconditioned gradient that filters out unstable fluctuations (Ke, You & Wang Reference Ke, You and Wang2026). While the eddy-viscosity approach does not reproduce the detailed structure of individual adjoint trajectories, it captures the large-scale sensitivity structure, which may be useful for interpreting the dependence of measurements.
A comparison in spectral space is presented in figure 6. The RA-ANS prediction (a) and the ensemble-averaged adjoint spectra (b) exhibit similar large-scale organisation, with energy concentrated in comparable regions of
$(k_x^+,k_z^+)$
space. At long times, the low wavenumbers show good agreement. These modes correspond to the dominant large-scale structures that are most robust under ensemble averaging. In contrast, the high-wavenumber components, which are difficult to converge in the ensemble-averaged adjoint approach, show significant discrepancies in their energy. This discrepancy reflects the strong sensitivity of small-scale adjoint structures to chaos, which also leads to slow statistical convergence.
Energy spectra of the DOD as a function of
$\boldsymbol{k}^+ = (k_x^+, k_z^+)$
, integrated in the wall-normal direction for
$Re_{\tau }=180$
. The energy is normalised by the maximum value among all wavenumber pairs. (a) normalised energy spectra from RA-ANS. (b) normalised energy spectra from averaging
$5 \times 10^3$
samples of ANS equations.
(a) Convergence of the ensemble adjoint method. Squared norm of the ensemble-averaged adjoint fields (
$\|\langle \boldsymbol{u}^{\dagger } \rangle \|^2$
, red solid lines), and ensemble average of the adjoint energy (
$\langle ||\boldsymbol{u}^{\dagger }||^2 \rangle$
, black solid lines). Thicker lines indicate
$\{1, 2, 3, 4, 5\} \times 10^3$
samples. Dashed black line marks the best exponential fit
$e^{0.0134\zeta ^+}$
for the averaged adjoint energy. The mean energy decay using RA-ANS is marked by the red dash-dotted line. (b) Ratio between the energy of the mean and mean energy for the adjoint field using
$5 \times 10^3$
samples.
Figure 7. Long description
Panel A: A line graph shows the convergence of the ensemble adjoint method. The x-axis is labeled with the variable zeta plus and ranges from 0 to 50. The y-axis is labeled with the squared norm of the ensemble-averaged adjoint fields and the ensemble average of the adjoint energy, ranging from 10^-8 to 10^-4. The graph includes red solid lines representing the squared norm of the ensemble-averaged adjoint fields, black solid lines representing the ensemble average of the adjoint energy, and a dashed black line marking the best exponential fit for the averaged adjoint energy. Thicker lines indicate larger sample sizes. A red dash-dotted line marks the mean energy decay using RA-ANS. Panel B: A line graph shows the ratio between the energy of the mean and mean energy for the adjoint field using samples. The x-axis is labeled with the variable zeta plus and ranges from 0 to 50. The y-axis is labeled with the ratio of the energy of the mean to the mean energy, ranging from 10^-4 to 10^0. The graph shows a decreasing trend.
In figure 7(
$a$
), we report the norm of the ensemble-averaged adjoint field (red) and of their energy (black), plotted versus backward time. Convergence is shown with the number of realisations (thin-to-thick lines). Even though each adjoint trajectory is amplifying exponentially, the ensemble-averaged adjoint field from five thousand realisations decays through
$\zeta ^+ \approx 50$
. A smaller number of samples is insufficient to yield converged statistics throughout the reported backward-time horizon. The RA-ANS prediction of the energy
$\|\overline {u}^{\dagger }\|^2$
is shown by the red dash-dotted line, and exhibits convincing agreement with the converged ensemble-averaged adjoint energy
$\|\left \langle \boldsymbol{u}^\dagger \right \rangle \|^2$
over the time interval
$20 \le \zeta ^+ \le 50$
. This agreement indicates that the model captures the dominant large-scale sensitivity despite the chaotic growth of individual realisations. In contrast, the ensemble-averaged energy
$\langle \|{\boldsymbol{u}^\dagger }\|^2\rangle$
grows exponentially with time, at twice the Lyapunov exponent. In panel (
$b$
), we plot the ratio of the two quantities,
$\|\left \langle \boldsymbol{u}^{\dagger } \right \rangle \|^2 / \left \langle ||\boldsymbol{u}^{\dagger }||^2 \right \rangle$
, which shows an exponential decay, thus indicating that convergence of the ensemble average becomes exponentially more difficult for longer backward times.
We verified that the RA-ANS predictions are not strongly influenced by the magnitude of the eddy viscosity by performing additional computations in which
$\nu _t(y)$
was uniformly scaled by factors of
$0.8$
and
$1.2$
relative to the nominal Cess profile. The resulting mean DOD fields exhibit no noticeable qualitative differences from the baseline case. This test supports the view that the dominant spatio-temporal organisation of the mean DOD is primarily governed by the underlying mean-shear dynamics, rather than by the precise magnitude of the turbulent diffusivity.
Forward evolution of DOD of a late-time streamwise wall-shear measurement (
$ t_m^+ = 200$
) in a turbulent channel with
$Re_{\tau } = 1000$
, shown using isosurfaces of streamwise velocity perturbation. Measurement location is marked by a vertical dashed line.
(a) The normalised energy spectra at three selected times
$t^+=\{40,120,200\}$
. Regions
$A$
and
$B$
denote the spectral ranges associated with the Orr mechanism and streaky structures, respectively. (b) Forward evolution of DOD of a late-time measurement (
$ t_m^+ = 200$
). Contours of (b.i) spanwise and (b.ii) streamwise averages of
$\overline {u}$
, illustrating the contribution of the Orr mechanism and streaky structures to the perturbation dynamics in a turbulent channel of
$Re_{\tau } = 1000$
. White dashed lines in the side view mark the inclination angle of the structure.
Figure 9. Long description
Panel A: Three heat maps show the normalized energy spectra at selected times. The x-axis represents the streamwise wavenumber kx+ and the y-axis represents the spanwise wavenumber kz+. The color scale indicates the logarithm of the normalized energy. Regions A and B denote the spectral ranges associated with the Orr mechanism and streaky structures, respectively. Panel B: A series of contour plots illustrate the forward evolution of the direct optimal detection of a late-time measurement. Subpanel (bi) shows the spanwise average of the perturbation dynamics in a turbulent channel, with the x-axis representing the streamwise direction x and the y-axis representing the wall-normal direction y. Subpanel (bii) shows the streamwise average, with the x-axis representing the wall-normal direction z. White dashed lines in the side view mark the inclination angle of the structure.
3.2. Forward dynamics of the mean DOD
To understand the dynamics of the mean DOD of a wall-stress measurement, we initialise a forward simulation with the adjoint state from a measurement at time
$t_m^+ = 200$
, by setting
$\overline {\boldsymbol{u}}_0 = \overline {\boldsymbol{u}}^\dagger$
. We then compute the forward evolution using the RA-LNS, and the results are reported in figure 8. The DOD can be clearly separated into two components: the first comprises upstream-tilted structures and the second features streaky structures. The upstream-tilted structures evolve rapidly under the influence of the mean shear, and the dependence of the wall-stress measurement on Orr amplification in the earlier flow. The corresponding wavenumbers in the energy spectra are identified as region
$A$
in figure 9(
$a$
). The peak streamwise wavenumber
$k_x^+$
remains nearly constant over time, consistent with localised non-modal amplification processes. On the other hand, the streaky structures have small spanwise wavenumber
$k_z^+$
, and are identified in the spectra (figure 9
$a$
) as region
$B$
. The peak
$k_z^+$
becomes smaller throughout the forward evolution, indicating an elongation of the streaks. Such structures align with the most energetic modes of the forward dynamics. Together, these two parts combine in an optimal manner to influence the measurement data at the measurement time.
To further distinguish these two components, we isolate each of them by averaging the evolved fields in the spanwise figure 9(
$b.i$
) and streamwise figure 9(
$b.ii$
) directions. This averaging effectively spectrally filters the field along the
$k_z^+=0$
and
$k_x^+=0$
axes, respectively. These two cases highlight the dominant contributions associated with the Orr-type structures and streaky motions. The spanwise-averaged
$ [ \overline {u} ]_z$
in figure 9
$(b.i)$
exhibits initial alignment against the mean shear and gradually rotates towards vertical alignment as time evolves, which is consistent with the Orr mechanism. The streamwise-averaged
$ [ \overline {u} ]_x$
in figure 9
$(b.ii)$
shows a self-similar expansion while maintaining its overall shape, characteristic of a streak growth dynamics.
For the Orr amplification, we compute the dominant orientation of the flow field from the phase of the streamwise velocity in Fourier space. Specifically, we consider the complex-valued mode shape
$\hat {\overline {u}}(y,t)$
associated with the wavenumber
$K^+_x$
corresponding to the maximum energy along the line
$k^+_z=0$
. The local phase is then defined as
$\phi (y,t) = \arg (\hat {\overline {u}}(y,t) )$
. The wall-normal wavenumber is obtained from the phase gradient as
$K^+_y = \partial \phi / \partial y$
, which together with
$K^+_x$
defines the local wavevector
$\boldsymbol{K}^+ = (K^+_x, K^+_y)$
. The inclination angle is defined as
\begin{equation} \theta =\frac { \int _0^{y_c}\arctan \left (\frac {K^+_y}{K^+_x}\right )||\hat {\overline {u}}(y,t)||^2 {\rm d}y}{\int _0^{y_c} ||\hat {\overline {u}}(y,t)||^2 {\rm d}y}, \end{equation}
using
$||\hat {\overline {u}}||^2$
as weights to suppress contributions from low-amplitude regions, with integration from the wall to the cutoff range
$y_c = 0.2$
. This procedure yields a time-resolved estimate of the dominant orientation of
$ [ \overline {u} ]_z$
in the
$(x,y)$
plane. As shown in figure 10
$(a)$
, the extracted orientation angles from forward evolution of DOD from two different measurement times exhibit a similar temporal dynamics: the structure is initially upstream tilted relative to the mean shear, then progressively rotates toward a near-vertical alignment at the measurement time. The time dependence of the wall-shear stress at the measurement locations,
$\overline {m} = \nu ({\partial \overline {u}(t)}/{\partial y} ) \vert _{\boldsymbol{x}_m}$
, is plotted in figure 10(b) normalised by the values at the measurement times,
$t_m^+ = 200$
and
$t_m^+ = 160$
. The curves reflect the transient amplification of the fields as they approach the measurement times, where the peak stress is recorded. The earlier-in-time amplification takes place as the adjoint structure undergoes reorientation by the mean shear, consistent with the Orr mechanism. These observations are qualitatively similar to the findings by Encinar & Jiménez (Reference Encinar and Jiménez2020), where conditionally averaged ejection and sweep events display a rotation from upstream-tilt to near-vertical configurations. Note, however, that the adjoint DOD computes the optimal precursor structure to maximise the measurement data (Wang & Zaki Reference Wang and Zaki2025b
), here, the wall-shear stress and not the perturbation kinetic energy. As such, there is no a priori justification for an agreement. For example, unlike energy-based optimal growth, the peak amplification occurs near
$\theta \approx -30^\circ$
, shortly after the structure passes through vertical alignment. The difference is amply justified by the difference in the objectives, and arguably more important is the persistent importance of Orr amplification.
(a) Temporal evolution of the near-wall orientation angle
$\theta$
, computed from the local phase gradient of
$[\bar {u}]_z$
. (b) Wall-shear stress at the sensor location, normalised by its value at
$t_m$
, from forward evolution of the DOD associated with (solid)
$t_m^+ = 200$
and (dashed)
$t_m^+ = 160$
.
The two-component decomposition of the mean adjoint field is effectively a distinction of two types of information contained in wall-shear-stress measurements, each associated with a different physical mechanism. Together, Orr amplification and lift-up define the flow regions and mechanisms, or information, that are observable from the measurements in an averaged sense, which are typically obscure in the chaotic dynamics of individual realisations of the adjoint fields.
4. Conclusion
Estimating the DOD of a measurement in turbulent channel flow is challenging due to the exponential growth of adjoint fields in backward time, which reflects the inherent difficulty of data assimilation in turbulence. While ensemble averaging of adjoint trajectories can, in theory, yield the mean DOD, this approach suffers from prohibitively slow convergence and overwhelming computational cost. Specifically, while the ensemble-averaged adjoint field should decay, the energy of individual adjoint realisations amplifies at twice the Lyapunov exponent and, as such, a very large number of adjoint trajectories are required in the ensemble averaging for convergence. We overcome these obstacles by introducing a linear eddy-viscosity closure to approximate the forward–adjoint correlation in the RA-ANS, thus enabling a direct and efficient computation of the mean DOD. This model bypasses the exponential divergence of individual adjoint trajectories, and its accuracy was demonstrated against the expensive ensemble average sensitivity up to a backward time
$\zeta ^+=\mathcal{O}(40)$
.
The present work establishes a connection between several concepts that are typically considered separately: adjoint-based sensitivity analysis, the Lyapunov-driven divergence of chaotic trajectories and the non-normal transient growth mechanisms that organise turbulent structures. In particular, the results highlight how the loss of information associated with long adjoint time horizons coexists with a coherent ensemble-averaged sensitivity that can be captured through appropriate regularisation using eddy viscosity for the ensemble-averaged adjoint equations.
Our results reveal that the mean DOD exhibits a universal structure across a wide range of Reynolds numbers, which comprises two distinct dynamical components: an upstream-tilted part associated with the Orr mechanism, which is rapidly reoriented by the mean shear to influence wall observations, and streamwise-elongated streaks that reflect optimal energy growth over longer time scales. These components combine in a dynamically optimal way to shape the DOD.
The DOD displays a marked asymmetry relative to its forward counterpart, the DOI, reflecting the non-self-adjoint nature of the governing operator in the presence of Reynolds stresses. This asymmetry highlights the limitations of interpreting observations based solely on the forward dynamics. In contrast, the mean DOD provides a more appropriate framework for identifying the flow structures that influence a given measurement, as directly quantifies the backward-time sensitivity to earlier flow states. Furthermore, the spatio-temporal and spectral characteristics of the DOD are important because they describe the information content in wall-stress measurement, which can be extracted using data assimilation.
Funding
The authors acknowledge financial support from the National Science Foundation (grant 2332057) and the Office of Naval Research (grant N00014-25-1-2290).
Declaration of interests
The authors report no conflict of interest.
Reτ=103
∂u/∂y
y=0
(kx+,kz+)
u¯†/u¯max†=±0.01
⟨u†⟩/⟨u†⟩max=±0.01
ζ+={4,8,12,16,20,50}
Reτ=180
u†/umax†=±0.01
ζ+={4,8,12,16,20,50}
Reτ=180
u¯†
⟨u†⟩
5000
u(n)†
0.7
k+=(kx+,kz+)
Reτ=180
5×103
‖⟨u†⟩‖2
⟨||u†||2⟩
{1,2,3,4,5}×103
e0.0134ζ+
5×103
tm+=200
Reτ=1000
t+={40,120,200}
A
B
tm+=200
u¯
Reτ=1000
θ
[u¯]z
tm
tm+=200
tm+=160
Open access
