1. Introduction
In turbulent channel flow, the walls play an important role in the generation of vorticity. Despite this role, whether wall observations of the stresses can be used in data assimilation to accurately predict the initial state of the flow and its evolution are challenging questions. The difficulty is primarily due to the nonlinear, chaotic nature of turbulence: small deviations in the initial state grow exponentially in time (Deissler Reference Deissler1986) and, conversely, small errors in the observations can obscure the reconstruction of the flow. Previous studies have demonstrated that the near-wall turbulence and only the outer large-scale structures can be reconstructed from wall signals. In this work, we adopt an adjoint-variational approach to quantify the domain of dependence of wall sensors, and demystify the difficulty of flow reconstruction from wall observations.
Previous attempts to estimate the state of wall-bounded turbulence from limited observations can be separated into two general classes: filters and smoothers (e.g. Colburn, Cessna & Bewley Reference Colburn, Cessna and Bewley2011; Suzuki Reference Suzuki2012; Mons et al. Reference Mons, Chassaing, Gomez and Sagaut2016; Wang, Hasegawa & Zaki Reference Wang, Hasegawa and Zaki2019b). In filtering techniques, the state is adjusted to best match the most recently available observations. On the other hand, smoothers take into account the entire observation history.
Extended (Hœpffner et al. Reference Hœpffner, Chevalier, Bewley and Henningson2005) and ensemble (Chevalier et al. Reference Chevalier, Hœpffner, Bewley and Henningson2006) Kalman filters were adopted to estimate turbulent channel flow from wall observations, and in both cases accuracy deteriorates with distance from the wall. Suzuki & Hasegawa (Reference Suzuki and Hasegawa2017) adopted linear stochastic estimation (LSE) to interpret wall observations at the friction Reynolds number of 
$Re_\tau = 100$. Accuracy of their predicted flow fields was commensurate with earlier efforts that used Kalman filtering techniques: the estimation was only successful in the region 
$y^{+} \lesssim 20$, where 
$+$ indicates scaling in viscous wall units, and the error in the estimated field increased beyond that height. At higher friction Reynolds number (
$Re_\tau = 1000\text {--}5000$), LSE was able to also capture the large-scale structures in the outer layer (Encinar & Jiménez Reference Encinar and Jiménez2019). It is important to note that LSE relies on prior knowledge of the correlation between wall observations and the flow, which may not be available. In addition, LSE and also Kalman filtering techniques do not satisfy the Navier–Stokes equations.
In contrast to filters, smoothers satisfy the governing equations and attempt to reproduce the observations over the entire time horizon during the forward evolution of the flow. The problem is formulated as a nonlinear optimization: an initial state is sought to minimize a cost function that is proportional to the difference between predicted and available observations. Both ensemble-variational (Mons, Wang & Zaki Reference Mons, Wang and Zaki2019; Buchta & Zaki Reference Buchta and Zaki2021) and four-dimensional adjoint-variational (4DVar) (Le Dimet & Talagrand Reference Le Dimet and Talagrand1986; Li et al. Reference Li, Zhang, Dong and Abdullah2020) data assimilation can be adopted for the minimization procedure. In the latter approach, which will be adopted herein, deviations from the observations appear in the adjoint system as a forcing term; the adjoint variable at the initial time provides the gradient of the cost function with respect to the initial state, and is used to adjust the initial condition. Adjoint-variational methods are popular in weather prediction (Kleist & Ide Reference Kleist and Ide2015a,Reference Kleist and Ideb), nonlinear stability analysis (Schmid Reference Schmid2007; Kerswell Reference Kerswell2018) and optimal flow control (Luchini & Bottaro Reference Luchini and Bottaro2014). In wall turbulence, Bewley & Protas (Reference Bewley and Protas2004) used 4DVar to reconstruct the initial state of channel flow from observations of wall friction and pressure, at a moderate friction Reynolds number 
$Re_{\tau } = 100$. The estimation was accurate near the wall, but was nearly uncorrelated with the true initial flow state in the channel centre. As the Reynolds number is increased, the outer large-scale structures can also be reconstructed from wall observations using 4DVar (Wang & Zaki Reference Wang and Zaki2021).
A measure of the difficulty of reconstructing the turbulent state from wall data is desirable. In the present work, we will objectively and quantitatively relate the difficulty of variational flow estimation to the sensitivity of wall measurements to the flow state. In variational methods, since the optimal initial condition minimizes the cost function, the local gradient vanishes. Therefore, the second-order sensitivity of the cost function with respect to the initial condition, or its Hessian matrix, characterizes the local geometric properties of the optimization problem. Specifically, the eigenvalues of the Hessian are the curvatures of the cost function if the initial condition were to be disturbed along the directions given by the associated eigenvectors. Higher eigenvalues indicate higher curvature, and therefore the data assimilation algorithm will predict these eigendirections from the measurements with higher certainty. The flow structures associated with smaller eigenvalues, or curvatures, are more difficult to estimate accurately. In addition, the ratio of the largest to the smallest eigenvalue of the Hessian is its condition number, which provides a quantitative measure of the ill-conditioned nature of the inverse problem.
The Hessian analysis has broad applications and profound implications. For example, Kalmikov & Heimbach (Reference Kalmikov and Heimbach2014) estimated the large-scale ocean state using an adjoint method, evaluated the Hessian matrix using algorithmic differentiation, and used its leading eigenvectors to study the influence of uncertainties in observations on estimation of oceanographic target quantities. In data assimilation, the landscape of the cost function depends on the design of measurements (e.g. choice of observed quantity, sensor placement, relative weighting of observations), and the Hessian analysis can provide an objective and quantitative approach for optimizing these measurements parameters in order to mitigate the difficulty of state estimation (Zaki & Wang Reference Zaki and Wang2021). For example, Mons et al. (Reference Mons, Wang and Zaki2019) optimized the placement of sensors for estimating scalar sources in turbulent channel flow by minimizing the condition number of the Hessian, and Buchta & Zaki (Reference Buchta and Zaki2021) adopted a Hessian-based strategy to optimize the weighting of wall-pressure measurements for estimating hypersonic transitional boundary-layer flow. The primary focus of the present work is on the role of the Hessian matrix in adjoint-variational estimation of turbulence from wall observations. In this context, the Hessian matrix has never been evaluated or analysed before – a gap that we address herein.
The adjoint-variational approach for state estimation is reviewed briefly in § 2.1, and sample flow reconstructions in turbulent channel flow are presented. Using the forward-adjoint duality relation, we formulate an approach for computing the Hessian matrix at the true initial flow state in § 2.2. The adjoint fields associated with wall observations are reported in § 3.1, and interpreted as the dependence of wall data on the flow at different times from the observations. The adjoint fields are then used to construct the Hessian matrix. Eigenanalysis for the Hessian is presented in § 3.2, and used to explain the sensitivity of the wall data to various wall-normal locations in state estimation. We proceed to evaluate the sensitivity of wall observations to the most energetic flow structures that are obtained from a proper orthogonal decomposition (POD) of turbulent channel flow. Finally, statistical behaviour of the adjoint field at long times is discussed in § 3.3.
2. Methods
The flow configuration of interest is statistically stationary turbulent channel flow. The reference length 
$h$ is the channel half-height and the reference velocity is the bulk value 
$\mathcal {U}$. The flow is governed by the non-dimensional incompressible Navier–Stokes equations
\begin{gather} \frac{\partial \boldsymbol{U}}{\partial t} + \left(\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla}\right) \boldsymbol{U} ={-}\boldsymbol{\nabla} P + \frac{1}{Re}\nabla^{2} \boldsymbol{U}, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{U} = 0, \end{gather}
\begin{gather}\left.\boldsymbol{U}\right\rvert_{t=0} = \boldsymbol{U}_0, \end{gather}
where 
$Re \equiv \mathcal {U} h / \nu$ is the bulk Reynolds number and 
$\nu$ is the kinematic viscosity. The spatially and temporally dependent velocity and pressure are denoted 
$\boldsymbol {U}$ and 
$P$, and 
$\boldsymbol {U}_0$ is the initial state which is the target of reconstruction. The flow is driven by a constant streamwise pressure gradient, is periodic in the streamwise (
$x$) and spanwise (
$z$) directions, and satisfies the no-slip conditions at the bottom and top walls 
$y = \{0,2\}$. The governing equations (2.1) are solved using a fractional step method with a volume-flux formulation (Kim & Moin Reference Kim and Moin1985; Rosenfeld, Kawak & Vinokur Reference Rosenfeld, Kawak and Vinokur1991). The diffusion terms are discretized implicitly in time with Crank–Nicolson scheme while the nonlinear advection terms are treated explicitly using Adams–Bashforth. The pressure Poisson equation is solved by Fourier transforms in the periodic streamwise and spanwise directions, followed by a tridiagonal solver in the wall-normal direction. The numerical algorithm has been used extensively for direct numerical simulations of transitional (Zaki Reference Zaki2013) and turbulent flows (Lee, Sung & Zaki Reference Lee, Sung and Zaki2017).
The majority of the results are focused on Reynolds numbers 
$Re = 2800$ and 10 935; the corresponding friction Reynolds numbers are 
$Re_\tau \equiv u_\tau h / \nu =180$ and 
$Re_\tau =590$, where 
$u_\tau \equiv \sqrt {\nu ({\rm d}\bar {U}/{{\rm d}y})_{wall}}$ is the friction velocity evaluated from the mean wall shear stress, and an overbar denotes averaging in the homogeneous spatial directions and in time. We adopted a Cartesian grid with uniform spacing in both the streamwise and spanwise directions and hyperbolic stretching in the wall-normal coordinate. The dimensions of the computational domains and grid resolutions are reported in table 1, including for additional Reynolds numbers that are referenced in the text. The grid resolutions were selected based on the general criteria that 
$\Delta x^{+} \lesssim 10$, 
$\Delta z^{+} \lesssim 5$ and 
$\Delta y^{+}_{min} \lesssim 1$.
Table 1. Domain sizes and grid resolutions. Grid sizes are normalized by the viscous length scale, e.g. 
$\Delta x^{+} \equiv \Delta x\,u_{\tau } / \nu$.
Independent reference simulations were performed to generate the observation data that are then adopted in the state estimation procedure. The reference flow fields are therefore the hidden truth that will be used to quantify the accuracy of the estimated fields.
2.1. Adjoint-variational state estimation
Starting from wall observations only, we seek to evaluate an estimate of the true initial state of the flow 
$\boldsymbol {U}_0$, which will be denoted 
$\tilde {\boldsymbol {U}}_0$. Similar to 
$\boldsymbol {U}$, the estimated field 
$\tilde {\boldsymbol {U}}$ satisfies the Navier–Stokes equations (2.1) which are referred to as the forward model. Errors in 
$\tilde {\boldsymbol {U}}_0$ result in deviation of the associated model predictions from available observations, and the difference is used in the definition of our cost function that we aim to minimize. Throughout the present study, observations are only available at the top and bottom surfaces 
$S$ of the channel. The initial state is defined at 
$t=0$, and the final observation time is denoted as 
$t_m$. When observations are only available at one instant, 
$t_m$, the associated cost function for the streamwise wall stress is
\begin{equation} \mathcal{J}_u(\tilde{\boldsymbol{U}}_0;t_m) = \frac{1}{2S}\int_S \left( \frac{\partial\tilde{U}}{\partial y}-\frac{\partial U}{\partial y} \right)^{2}_{(\boldsymbol{x}_m,t_m)}\,{\rm d}x_m \,{\rm d}z_m, \end{equation}
where 
$\boldsymbol {x}_m$ is the observation location. Similarly, cost functions 
$\mathcal {J}_w$ and 
$\mathcal {J}_p$ are defined for deviations from instantaneous observations of spanwise wall stress 
$({\partial \tilde {W}}/{\partial y}-{\partial W}/{\partial y} )_{(\boldsymbol {x}_m,t_m)}$ and wall pressure 
$(\tilde {P}- P )_{(\boldsymbol {x}_m,t_m)}$.
Results for observations of only one component of the stress or pressure at one time instance are the focus of this work, and will be discussed in detail in §§ 3.1 and 3.2. In order to place that discussion in context, however, we start by examining the capacity for state estimation when all three components are available as a function of time during the assimilation window; this preliminary step is a summary of a recent detailed study by Wang & Zaki (Reference Wang and Zaki2021). The cost function for the collective observations over the entire time horizon is defined as
\begin{equation} J = \frac{1}{Re^{2}} J_u + \frac{1}{Re^{2}} J_w + J_p, \end{equation}
where 
$J_u$ is the integral of the instantaneous cost function (2.2) over the observation window,
\begin{equation} J_u(\tilde{\boldsymbol{U}}_0) = \int_0^{t_m} \mathcal{J}_u(\tilde{\boldsymbol{U}}_0;t^{\prime}) \,{\rm d}t^{\prime} = \frac{1}{2S}\int_0^{t_m}\int_S \left( \frac{\partial\tilde{U}}{\partial y}-\frac{\partial U}{\partial y} \right)^{2}_{(\boldsymbol{x}_m,t^{\prime})} \,{\rm d}x_m \,{\rm d}z_m \,{\rm d}t^{\prime}, \end{equation}
and similarly for 
$J_w$ and 
$J_p$. The coefficients of 
$J_u$ and 
$J_w$ in (2.3) ensure that the three terms are comparable in magnitude.
The gradient of the cost function with respect to the initial condition is obtained by solving the adjoint equations
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}^{*} ={-}\frac{\partial J}{\partial \tilde{P}}, \end{gather}
\begin{gather}\frac{\partial \boldsymbol{u}^{*}}{\partial \tau} + \left( \boldsymbol{\nabla} \tilde{\boldsymbol{U}} \right) \boldsymbol{\cdot} \boldsymbol{u}^{*} - \left(\tilde{\boldsymbol{U}} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u}^{*} = \boldsymbol{\nabla} p^{*} + \frac{1}{Re} \nabla^{2} \boldsymbol{u}^{*} + \frac{\partial J}{\partial \tilde{\boldsymbol{U}}}, \end{gather}
where 
$*$ denotes adjoint variables, and 
$\tau \equiv t_m-t$ is the reverse time over the duration of the assimilation window 
$t_m$. The boundary conditions for the adjoint variables are no-slip at the walls and periodicity in the horizontal 
$x$ and 
$z$ directions. The source terms on the right-hand side of the equations are due to the deviation between the estimated and true observations. The gradient of the cost function with respect to the initial condition is obtained at the end of the adjoint computation, 
${\partial J}/{\partial \tilde {\boldsymbol {U}}_0} = \boldsymbol {u}^{*}(\tau =t_m)$ (for a detailed derivation, see Wang, Wang & Zaki (Reference Wang, Wang and Zaki2019a)).
The adjoint equations (2.5) are solved using the discrete-adjoint approach, which provides a more accurate gradient of the cost function than the continuous-adjoint counterpart (Vishnampet, Bodony & Freund Reference Vishnampet, Bodony and Freund2015). The gradient is then used in a quasi-Newton method in order to update the estimation of the initial state; our particular choice of quasi-Newton method was the limited-memory Broyden–Fletcher–Goldfarb–Shanno (L-BFGS) algorithm (Nocedal Reference Nocedal1980), which does not require the Hessian of the cost function and instead approximates it from successive gradients (see e.g. chapter 6 in Nocedal & Wright (Reference Nocedal and Wright2006)). The entire forward-adjoint gradient-based optimization is then repeated until the estimated initial state significantly reduces the cost function.
In summary, the procedure for minimizing the cost function (2.2) is the following: (i) start with a guess of the initial condition 
$\tilde {\boldsymbol {U}}_0$, march the forward equations (2.1) from 
$t=0$ to 
$t=t_m$, and store the instantaneous fields 
$\tilde {\boldsymbol {U}}(t)$; (ii) solve the adjoint equations (2.5) from 
$\tau = 0$ to 
$\tau = t_m$ and compute the gradient of the cost function; (iii) update the estimated initial condition 
$\tilde {\boldsymbol {U}}_0$, and repeat procedures (i) and (ii) until convergence. Note that this optimization process is computationally expensive, both in terms of storage cost and computational time. The former is caused by the storage of the forward flow fields 
$\tilde {\boldsymbol {U}}(t)$, and the latter is due to the hundreds of forward and adjoint simulations that are required for the minimization of the cost function.
Due to the exponential divergence of state-space trajectories of both the forward and adjoint dynamical systems, the choice of the assimilation time horizon is guided by the Lyapunov time scale. In viscous units, that horizon is 
$t_m^{+} = 50$, which corresponds to approximately one Lyapunov time scale, 
$\tau _{\sigma }^{+} \approx 50$ at 
$Re_{\tau } = 180$, and fields were stored at every time step to perform adjoint simulations. At each Reynolds number, the initial guess of 
$\tilde {U}_0$ is obtained from LSE (Adrian & Moin Reference Adrian and Moin1988; Encinar & Jiménez Reference Encinar and Jiménez2019) using wall measurements. This guess is then refined using 100 L-BFGS iterations, each comprising a forward and an adjoint fully resolved simulations of turbulence over the time window of assimilation. The convergence history for different Reynolds numbers are similar, and the cost function is reduced to less than 
$4\,\%$ of its initial value.
A sample estimation of the turbulence from wall observations is provided in figure 1. Since the estimated flow progressively approaches the true state in time within the estimation window (e.g. Wang et al. Reference Wang, Wang and Zaki2019a), we only report the predictions at 
$t=t_m$ when the estimation is most accurate. Note that the temporal decay of the errors is unique to 4DVar. In contrast, starting from an LSE of the initial state and solving the Navier–Stokes equations would lead to a flow evolution that diverges exponentially from the true state-space trajectory. The contours in figure 1(a) show both the streamwise velocity fluctuations from the turbulence in the true state and in the reconstruction using our 4DVar algorithm, at 
$Re_{\tau } = 590$. The two fields are nearly identical very near the wall, but precipitously deviate from one another from the buffer layer and into the bulk of the channel. The turbulence within the core of the channel is not accurately reconstructed, with the exception of the large-scale streaky structures; these outer energetic motions are known to modulate the near-wall region and the wall shear stress (Abe, Kawamura & Choi Reference Abe, Kawamura and Choi2004; Mathis, Hutchins & Marusic Reference Mathis, Hutchins and Marusic2009; Hwang et al. Reference Hwang, Lee, Sung and Zaki2016; You & Zaki Reference You and Zaki2019). The ability to predict the outer large scales from wall data using 4DVar is qualitatively consistent with, but quantitatively more accurate than, other estimation techniques, for example applying LSE at 
$t=t_m$ (Suzuki & Hasegawa Reference Suzuki and Hasegawa2017; Encinar & Jiménez Reference Encinar and Jiménez2019). The herein reported characteristics of the estimated field using 4DVar are a summary of the findings by Wang & Zaki (Reference Wang and Zaki2021).
Figure 1. (a) Top views of streamwise velocity fluctuations in the true and estimated states at 
$t=t_m$ for 
$Re_{\tau } = 590$ at selected 
$y$ locations, 
$y^{+} = \{15, 100, 300\}$. (b) Correlation coefficient between true and estimated streamwise velocity fluctuations at 
$t=t_m$. Progressively darker lines corresponding to 
$Re_{\tau } = \{100,180,392\}$; thick blue curve corresponds to 
$Re_{\tau } = 590$; filled circles mark the wall-normal locations of the subpanels in panel (a).
The correlation coefficient 
$\mathcal {C}_{xz}$ between the true and estimated fields at the final time, 
$t=t_m$, is reported in figure 1(b),
\begin{equation} \mathcal{C}_{xz} = \frac{ \langle \tilde{U}^{\prime} U^{\prime} \rangle_{xz} } { { \langle \tilde{U}^{\prime 2} \rangle_{xz}^{1/2}} {\langle U^{\prime 2} \rangle_{xz}^{1/2}} }, \end{equation}
where angle brackets denote averaging with respect to the marked dimensions and prime indicates fluctuations 
$U^{\prime } = U - \langle U\rangle _{xz}$. The correlation coefficient decays from the start of the buffer layer and the estimated field becomes essentially uncorrelated with the true state as we approach the channel centre. We have repeated the analysis for Reynolds numbers 
$Re_\tau = \{100, 180, 392, 590\}$, and the findings remain qualitatively similar. In effect, as the Reynolds number increases, wall observations can be interpreted to accurately estimate the flow in a diminishingly small physical region near the wall and also the outer large-scale structures that are known to have a near-wall signature.
2.2. Hessian analysis
A fundamental difficulty in reconstructing the initial condition is due to the chaotic nature of the governing equations and their adjoint: just as infinitesimal deviation in the initial conditions leads to exponentially diverging state-space trajectories in forward time, small mismatch between the estimated and true observations leads to exponentially diverging trajectories in the adjoint. In terms of the variational state-estimation procedure, we can quantify this difficulty by examining the convergence behaviour when the estimated state 
$\tilde {\boldsymbol {U}}_0$ is infinitesimally close to the true solution 
$\boldsymbol {U}_0$. We can then introduce the deviation field 
$\boldsymbol {u} = \tilde {\boldsymbol {U}}-\boldsymbol {U}$, which is governed by the linearized Navier–Stokes equations,
\begin{gather} \frac{\partial \boldsymbol{u}}{\partial t} + \left(\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u} + \left(\boldsymbol{u} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{U} ={-} \boldsymbol{\nabla} p + \frac{1}{Re} \nabla^{2} \boldsymbol{u}, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = \boldsymbol{0}. \end{gather}
The initial condition of the linearized equations is 
$\boldsymbol {u}_0 = \tilde {\boldsymbol {U}}_0 - \boldsymbol {U}_0$.
Here, we focus on the estimation of the initial condition using the instantaneous cost function (2.2). In other words, we only consider instantaneous observations at 
$t=t_m$, recorded over the entire walls, and in this manner we can evaluate the influence of the observation time 
$t_m$. Extension to the full assimilation window involves integration over time, which is straightforward. In terms of the new variable 
$u$, the cost function (2.2) becomes
\begin{equation} \mathcal{J}_u(t_m) = \frac{1}{2S}\int_S \left( \frac{\partial u}{\partial y} \right)^{2}_{(\boldsymbol{x}_m,t_m)} \,{\rm d}x_m \,{\rm d}z_m, \end{equation}
and similar expressions can be written for 
$\mathcal {J}_w$ and 
$\mathcal {J}_p$. Note that at the true solution 
$\boldsymbol {u}_0 = \tilde {\boldsymbol {U}}_0 - \boldsymbol {U}_0=\boldsymbol {0}$, and the gradient of the cost function (generically denoted as 
$\mathcal {J}$) vanishes,
\begin{equation} \left.\frac{\partial \mathcal{J}}{\partial \tilde{\boldsymbol{U}}_0}\right\rvert_{\tilde{\boldsymbol{U}}_0 = \boldsymbol{U}_0} = \left.\frac{\partial \mathcal{J}}{\partial \boldsymbol{u}_0} \right\rvert_{\boldsymbol{u}_0 = \boldsymbol{0}} = \boldsymbol{0}. \end{equation}Therefore, the difficulty of the state estimation is directly tied to the second-order derivative, or Hessian matrix, of the cost function. If the Hessian is well conditioned, converging from the vicinity of the true solution will be straightforward; conversely, an ill-conditioned Hessian obstructs the optimization algorithm from finding the true state, even if the initial guess is infinitesimally close to the true solution.
The cost function 
$\mathcal {J}$ depends on 
$\boldsymbol {u}$, and involves solving the linearized forward Navier–Stokes equations using the initial condition 
$\boldsymbol {u}_0$ and registering observations. Computing the Hessian using this forward relation is very expensive, since each possible spatial location in the initial condition and each velocity component must be perturbed independently, and therefore 
$O(N_x \times N_y \times N_z \times 3)$ simulations are required, each evolved up to the time of the observation 
$t_m$. A more efficient method is sought using the duality relation of the linearized Navier–Stokes equations and their adjoint.
Our proposed approach for evaluating the Hessian matrix starts from the impact of a deviation in the velocity field, 
$\boldsymbol {u}$, on an isolated observation at 
$(\boldsymbol {x}_m,t_m)$. Mathematically, this quantity can be expressed as 
$[\boldsymbol {u}(t = t_m), \boldsymbol {\phi }(\boldsymbol {x}_m)]$, where 
$[\boldsymbol {a},\boldsymbol {b}] =\int _V \boldsymbol {a}^{T}\boldsymbol {b} \,{\rm d}V$ is the spatial inner product evaluated over the entire channel volume, and 
$\boldsymbol {\phi }$ is the observation kernel that extracts the measurements of interest. For example, when observing the streamwise wall-shear-stress, the kernel 
$\boldsymbol {\phi }$ is the wall-normal derivative, selectively acting on the streamwise component of 
$\boldsymbol {u}$ and evaluated at 
$\boldsymbol {x}_m$. Note that observing the deviation in wall pressure can also be expressed as 
$[\boldsymbol {u}(t = t_m), \boldsymbol {\phi }(\boldsymbol {x}_m)]$ because 
$p$ is linearly related to 
$\boldsymbol {u}$ by the linearized Navier–Stokes equations (2.7). Now we exploit forward-adjoint duality to rewrite this expression in terms of the adjoint field 
$\boldsymbol {u}^{*}$, which will simplify the evaluation of the Hessian,
\begin{align} {[}\boldsymbol{u}(t = t_m), \boldsymbol{\phi}(\boldsymbol{x}_m)] = [\mathcal{A}\boldsymbol{u}_0, \boldsymbol{\phi}(\boldsymbol{x}_m)] = [\boldsymbol{u}_0, \mathcal{A}^{*}\boldsymbol{\phi}(\boldsymbol{x}_m)] = [\boldsymbol{u}_0, \boldsymbol{u}^{*}(\tau = t_m;\boldsymbol{x}_m,t_m)]. \end{align}
In the first equality, we replaced 
$\boldsymbol {u} = \mathcal {A} \boldsymbol {u}_0$, where 
$\mathcal {A}$ is the forward linearized Navier–Stokes operator (2.7) that advances the initial condition 
$\boldsymbol {u}_0$ to the final state 
$\boldsymbol {u}$ at time 
$t = t_m$. The second equality introduces the adjoint operator 
$\mathcal {A}^{*}$; its action 
$\mathcal {A}^{*}\boldsymbol {\phi }$ is equivalent to solving the linearized adjoint equations
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u}^{*} = 0, \end{gather}
\begin{gather}\frac{\partial \boldsymbol{u}^{*}}{\partial \tau} +\left(\boldsymbol{\nabla} \boldsymbol{U}\right) \boldsymbol{\cdot} \boldsymbol{u}^{*} - \left(\boldsymbol{U} \boldsymbol{\cdot} \boldsymbol{\nabla} \right) \boldsymbol{u}^{*} = \boldsymbol{\nabla} p^{*} + \frac{1}{Re} \nabla^{2} \boldsymbol{u}^{*}, \end{gather}
\begin{gather}\boldsymbol{u}^{*}(\tau = 0) = \boldsymbol{\phi}(\boldsymbol{x}_m), \end{gather}
where 
$\tau = t_m-t$ is the reverse time, and the interval 
$0\le t \le t_m$ is the same time horizon as the forward equations. The adjoint boundary conditions are, as before, no-slip at the walls and periodicity in the horizontal directions. It is important, however, to note the difference between these adjoint equations (2.11) and the earlier ones (2.5); here the equations do not feature a forcing term by the cost function and their initial condition is the observation kernel. The resulting adjoint field is therefore not the gradient of the cost function, and its meaning will be explained below in terms of the duality relation (2.10). In the last equality of (2.10), we replace 
$\mathcal {A}^{*}\boldsymbol {\phi }$ by the resulting field at the end of the adjoint evolution, 
$\boldsymbol {u}^{*}(\tau =t_m)$.
Before continuing the derivation of the Hessian, we note that the duality relation (2.10) is a powerful expression: instead of evolving different initial deviations 
$\boldsymbol {u}_0$ and evaluating their respective impact on the observations, we can perform a single adjoint simulation starting from the observation kernel 
$\boldsymbol {\phi }$ and determine the impact on the observations by the simple inner product, 
$[\boldsymbol {u}_0, \boldsymbol {u}^{*}(\tau = t_m;\boldsymbol {x}_m,t_m)]$, for any choice of 
$\boldsymbol {u}_0$ and without performing any forward simulations. From this expression, if 
$\boldsymbol {u}^{*}$ is non-zero only within a compact region, the initial perturbation must be finite within the same region in order to affect observations. In other words, the support of 
$\boldsymbol {u}^{*}$ represents the domain of dependence of the observation evaluated at 
$(\boldsymbol {x}_m, t_m)$.
We now complete the derivation of the Hessian. Using the duality relation (2.10), we can rewrite the cost function with observation kernel 
$\boldsymbol {\phi }$ at time 
$t_m$ in terms of adjoint field 
$\boldsymbol {u}^{*}$,
\begin{equation} \mathcal{J}(\boldsymbol{u}_0;t_m) = \frac{1}{2S}\int_S [\boldsymbol{u}, \boldsymbol{\phi}]^{2} {\rm d}x_m \,{\rm d}z_m = \frac{1}{2S}\int_S [\boldsymbol{u}_0, \boldsymbol{u}^{*}]^{2} {\rm d}x_m \,{\rm d}z_m. \end{equation}
The quadratic relation between 
$\mathcal {J}$ and 
$\boldsymbol {u}_0$ thus becomes explicit, and facilitates the evaluation of both the gradient and the Hessian,
\begin{equation} \left.\begin{gathered} \frac{\partial \mathcal{J}}{\partial \boldsymbol{u}_0} = \frac{1}{S} \int_S\left[\boldsymbol{u}_0, \boldsymbol{u}^{*}\right] \boldsymbol{u}^{*} \,{\rm d}x_m \,{\rm d}z_m,\\ \boldsymbol{\mathcal{H}} \equiv \frac{\partial^{2} \mathcal{J}}{\partial \boldsymbol{u}_0\partial \boldsymbol{u}_0} = \frac{1}{S}\int_S \boldsymbol{u}^{*}\boldsymbol{u}^{*} \,{\rm d}x_m \,{\rm d}z_m. \end{gathered}\right\}\end{equation}
At optimality, 
$\boldsymbol {u}_0 = \boldsymbol {0}$ and the gradient vanishes while the Hessian matrix is finite. The integral in (2.13) is over all measurement locations, i.e. the top and bottom walls; the integrand is the cross-correlation of the adjoint fields, 
$\boldsymbol {u}^{*}\boldsymbol {u}^{*}$, which is the Hessian for an isolated measurement. In the language of control theory (Rowley Reference Rowley2005), the Hessian matrix corresponds to the observability Gramian. The size of 
$\boldsymbol {\mathcal {H}}$ is 
$[3N_xN_yN_z]^{2}$, where three is the number of velocity components and 
$N_xN_yN_z$ is the number of spatial grid points. Written explicitly,
\begin{equation} \mathcal{H}_{ij}(\boldsymbol{x}_1,\boldsymbol{x}_2; t_m) = \frac{1}{S}\int_S u_i^{*}(\boldsymbol{x}_1,\tau = t_m ;\boldsymbol{x}_m, t_m) u_j^{*}(\boldsymbol{x}_2, \tau = t_m;\boldsymbol{x}_m, t_m) \,{\rm d}x_m \,{\rm d}z_m. \end{equation}
From (2.14), and the fact that 
$\mathcal {J} \geq 0$, the Hessian at optimality must be symmetric and positive semidefinite. As a result, it has real eigenvalues and orthogonal eigenvectors. For each choice of the observation kernel 
$\boldsymbol {\phi }$, be it the wall-normal gradient or sampling of the local pressure, we define 
$\mathcal {J}_{u},\mathcal {J}_{w}$, 
$\mathcal {J}_{p}$ and use the above outlined procedure to evaluate the associated Hessian matrices, 
$\boldsymbol {\mathcal {H}}_u, \boldsymbol {\mathcal {H}}_w$, 
$\boldsymbol {\mathcal {H}}_p$.
The Hessian matrix in (2.14) is evaluated in physical space, has a very large size, and is difficult to compute and analyse. The integrand features the adjoint fields released from the measurement location 
$\boldsymbol {x}_m$ and evaluated at 
$\boldsymbol {x}_1$ and 
$\boldsymbol {x}_2$. Since we consider all measurement positions on the wall, and since the forward flow is homogeneous in the horizontal plane, the Hessian is invariant under a coordinate shift in 
$x$ and 
$z$ directions. Due to this homogeneity, the eigenvectors in 
$x$ and 
$z$ are guaranteed to be Fourier modes (similar to the eigenfunctions from POD, see e.g. Lumley (Reference Lumley1981)). As a result, after performing a Fourier transform of the adjoint fields, the Hessian can be compactly expressed in spectral space,
\begin{align} \hat{\boldsymbol{\mathcal{H}}}(k_x, k_z, y_1, y_2 ; t_m) =\int_S \hat{\boldsymbol{u}}^{*}(y_1,\tau = t_m;\boldsymbol{x}_m, t_m)\hat{\boldsymbol{u}}^{*}_{{\rm c.c}.}(y_2, \tau = t_m;\boldsymbol{x}_m, t_m) \,{\rm d}x_m \,{\rm d}z_m, \end{align}
where 
$\hat {\boldsymbol {u}}^{*}$ is the Fourier transform of 
$\boldsymbol {u}^{*}$ in 
$x$ and 
$z$, and subscript c.c. denotes a complex conjugate. For each wavenumber pair 
$(k_x, k_z)$, the Hessian has size 
$3N_y \times 3N_y$, and therefore its (eigen)analysis is affordable. Note that 
$\hat {\boldsymbol {\mathcal {H}}}$ preserves the same eigenvalues as the Hessian in (2.14), and the eigenfunctions of 
$\hat {\boldsymbol {\mathcal {H}}}$ at every 
$(k_x^{+},k_z^{+})$ are the spectral representation of those of the original Hessian (2.14). A major advantage of evaluating and analysing 
$\hat {\boldsymbol {\mathcal {H}}}$ is the relative computational efficiency.
In § 3.2, we will perform an eigenanalysis of 
$\hat {\boldsymbol {\mathcal {H}}}$ at each 
$(k_x,k_z)$ pair in order to extract the associated eigenmodes whose significance is threefold. Firstly, the eigenvalues are curvatures of the cost function along directions represented by the eigenvectors. Therefore, the performance of 4DVar state estimation can be explained in terms of the Hessian eigenmodes. Specifically, flow structures described by the leading modes will be most effectively targeted by the 4DVar procedure or, in other words, accurately reconstructed from observations. Secondly, the leading eigenvectors are also optimal perturbations to the forward field, that yield the largest influence on wall observations, and the corresponding sensitivity is quantified by the eigenvalues. Finally, we note that the condition number of the Hessian is a quantitative measure of the ill-posedness of the state-estimation problem.
The integral (2.15) is evaluated by sampling different measurement locations on the wall. It is also possible to use 
$\hat {\boldsymbol {u}}^{*}$ evaluated from different realizations of the true forward flow which is homogeneous in the horizontal dimensions and time. In total, 512 adjoint simulations were included in the evaluation of Hessian (2.15). This process was repeated for every observation time 
$t_m$ of interest and for every type of measurement: 
$\partial U/\partial y \vert _{wall}$, 
$\partial W/\partial y\vert _{wall}$ and 
$P\vert _{wall}$. The herein described approach for evaluating the Hessian matrix with the aid of the adjoint and in spectral space significantly reduces the computational cost relative to forward methods. Nonetheless, the adjoint simulations still require approximately 
$1.2 \times 10^{6}$ CPU hours for each Reynolds number considered.
3. Results
This section initially focuses on turbulent channel flows at 
$Re_{\tau } = 180$ and 
$590$. Ensemble averaged adjoint fields 
$\overline {\boldsymbol {u}^{*}}$ will be reported for different types of wall data and reverse times. Appropriate normalization demonstrates the similarity across different Reynolds numbers. The results also highlight the domains of dependence of different wall sensing modalities and their change with the time separation between the observation and the initial state. Eigendecomposition of the Hessian will be performed, and the eigenvalues and vectors at different wavenumbers and reverse times will provide a unique perspective on the problem of initial-state estimation. The section concludes with a discussion of the long-time statistical properties of the adjoint fields for select Reynolds numbers from table 1, 
$Re_{\tau } = \{180, 590, 1000\}$, and their implications.
3.1. Domain of dependence of isolated wall observation
The duality relation (2.10) provides an interpretation of the adjoint, where the support of 
$\boldsymbol {u}^{*}(\tau =t_m)$ can be viewed as the domain of dependence for an observation. Starting from an impulse at observation locations 
$\boldsymbol {x}_m$ on the wall, the adjoint field is a growing puff that propagates upstream. We evaluated the adjoint associated with different observation locations 
$\boldsymbol {x}_m$ on the wall, each evolving in backward time using the stored forward velocity fields (see (2.11)). For each observation location, we have also considered different starting times. Due to the homogeneity of channel-flow turbulence in the horizontal plane and time, the structures of all these adjoint fields 
$\boldsymbol {u}^{*}(\tau = t_m)$ are statistically similar within a coordinate shift. The ensemble average of the adjoint fields was evaluated and is reported in figure 2.
Figure 2. Three-dimensional isosurfaces of the ensemble-averaged, streamwise adjoint velocity fields for different measurement times 
$t_m^{+}$ when observing (a) 
${\partial U}/{\partial y}\vert _{wall}$, (b) 
${\partial W}/{\partial y}\vert _{wall}$ and (c) 
$P\vert _{wall}$, plotted at 
$\tau ^{+} = t_m^{+}$.
Three different types of observations are examined, namely 
$\partial U/\partial y \vert _{wall}$, 
$\partial W/\partial y\vert _{wall}$ and 
$P\vert _{wall}$. For each measurement modality, 512 samples were included in the ensemble average. Figure 2 shows that the fields are similar at the two Reynolds numbers, when visualized using viscous scaling. No matter the choice of wall observations, the adjoint structures all advect upstream of the observation point and expand as a function of reverse time, or for longer separation between the observation and the initial state. The quantity plotted is the streamwise adjoint velocity, specifically 
$u^{*}/u^{*}_{max}$. The spanwise symmetry when observing 
$\partial U/\partial y \vert _{wall}$ and 
$P\vert _{wall}$ can be contrasted to the antisymmetry when measuring 
$\partial W/\partial y\vert _{wall}$; in the first two cases, initial disturbances to the forward 
$U$ velocity at either side of the measurement cannot be distinguished, but they have opposite influences on a measurement of the spanwise stress. In addition, due to the non-local nature of pressure, the adjoint structure associated with wall-pressure observations, 
$P\vert _{wall}$, are appreciably larger than those associated with the wall stresses, especially for shorter times. This observation is congruent with the finding by Bewley & Protas (Reference Bewley and Protas2004) that including wall-pressure observations in addition to the surface stresses improves the accuracy of state estimation.
Figure 3 shows an instantaneous adjoint field for an even longer separation between the observation and the initial state, 
$\tau ^{+} = t_m^{+} = 220$ at 
$Re_{\tau } = 180$. The isosurfaces of adjoint streamwise velocity form a large chaotic patch that spans most of the horizontal plane, and which resembles the familiar turbulent spots from transitional flows (Cantwell, Coles & Dimotakis Reference Cantwell, Coles and Dimotakis1978; Marxen & Zaki Reference Marxen and Zaki2019) but with the opposite orientation and direction of propagation. The long-time behaviour of the adjoint field becomes independent of its initial condition, and is only meaningful to examine statistically. In addition, due to the chaotic nature of the underlying turbulent field 
$\boldsymbol {U}$, the convergence behaviour of the adjoint statistics quickly deteriorates and the required number of samples rapidly increases with larger 
$t_m$ (see e.g. Eyink, Haine & Lea Reference Eyink, Haine and Lea2004; Chandramoorthy et al. Reference Chandramoorthy, Fernandez, Talnikar and Wang2019). A discussion of the asymptotic behaviour of the adjoint field and its impact on the state estimation problem is provided in § 3.3, after we examine eigenproperties of the Hessian matrix at short times and the implication for state estimation (cf. figure 1).
Figure 3. Instantaneous, streamwise adjoint velocity field at 
$\tau ^{+} = t_m^{+} = 220$ due to an impulse of 
${\partial u}/{\partial y}\vert _{wall}$ at the wall for 
$Re_{\tau } = 180$. The isosurfaces show 
$u^{*} = \pm 0.001u_{max}^{*}$, and are coloured by the vertical distance from the wall.
3.2. Eigenvalues and eigenvectors of the Hessian
Recalling (2.15), the eigenmodes of the Hessian are Fourier modes in the homogeneous 
$x$ and 
$z$ directions, and their 
$y$-dependence can be evaluated from 
$\hat {\boldsymbol {\mathcal {H}}}$ for different wavenumber pairs 
$(k_x,k_z)$. With reference to the eigenvalues and vectors of the Hessian, we can then explain the capacity of 4DVar to reconstruct different wavenumber components of the initial flow state. The largest eigenvalues are the highest curvatures of the cost functions, and the associated eigenvectors are the flow structures which are most accurately estimated from the measurements.
The contours in figure 4 correspond to the largest eigenvalue of 
$\hat {\boldsymbol {\mathcal {H}}}_u$, the Hessian for instantaneous observation of 
$\partial U/\partial y \vert _{wall}$, as a function of the wavenumber vector. The results for 
$Re_{\tau } = 180$ and 
$590$ are unified by adopting viscous scaling, and a short observation time is considered, 
$t_m^{+} = 4$. The two Reynolds numbers 
$Re_{\tau } = 180$ and 
$590$ are compared using the colour and line contours, respectively, which demonstrates their agreement when scaled using viscous units. The maximum eigenvalue generally decays with increasing magnitude of the horizontal wavenumber vector 
$(k_x^{+}, k_z^{+})$, and the supremum corresponds to a two-dimensional mode 
$(k_x^{+}, k_z^{+}) = (0.12, 0)$. Figure 4(b) shows the entire eigenvalue spectra of 
$\hat {\boldsymbol {\mathcal {H}}}_u$ at that wavenumber pair, for both Reynolds numbers. The eigenvalues appear in pairs, and decay exponentially, which demonstrates that the sensitivity of the system is dominated by few directions associated with the leading eigenvectors: a typical feature of ill-conditioned systems.
Figure 4. (a) The logarithm of the largest eigenvalue of the Hessian matrix 
$\hat {\boldsymbol {\mathcal {H}}}_u$ at 
$t_m^{+} = 4$. Colour contours show the eigenvalues for 
$Re_{\tau } = 180$ and line contours correspond to 
$Re_{\tau } = 590$. (b) Eigenvalue spectra at the wavenumber pair 
$(k_x^{+},k_z^{+}) = (0.12, 0)$; the thicker line is the higher Reynolds number.
The profiles of select eigenvectors, 
$(\hat {u}(y), \hat {v}(y),\hat {w}(y))$, are examined in figure 5. The mode associated with the supremum eigenvalue at 
$(k_x^{+}, k_z^{+}) = (0.12,0)$ corresponds to spanwise homogeneous rolls. The roll peaks at 
$y^{+} \approx 3$ and generates counter vorticity at the wall that most effectively impacts the wall observation of streamwise shear stress. This eigenmode should not be interpreted using the conventional wisdom from forward simulations of turbulence and knowledge of its structures; the interpretation should be based on the notion of the adjoint. The eigenfunction shows the direction along which the cost function has the largest curvature, and hence the measurement is most sensitive to a perturbation in this direction at 
$\tau = t_m$ earlier in time. When this time is short, as is the case in figure 5, the forward dynamics do not have sufficient time to evolve an initial disturbance of the form of the eigenfunction into an amplified flow response that has a significant signature at the sensor. As such, the eigenfunction does not resemble flow structures that are common in simulations of turbulence. Rather the eigenfunction is the initial disturbance that, without significant development, essentially itself maximizes the sensor signal. It is nonetheless the leading direction in the process of minimizing the cost function during an adjoint-variational data assimilation using measurement of streamwise shear stress at 
$t_m^{+} = 4$.
Figure 5. Visualizations of the leading eigenmodes of the Hessian matrix with wavenumber 
$(k_x^{+},k_z^{+})=(0.12, 0)$ obtained for 
$t_m^{+} = 4$: (a) 
$Re_{\tau } = 180$; (b) 
$Re_{\tau } = 590$. Side view contours of the spanwise vorticity 
$\omega _z \equiv {\partial {v}}/{\partial x} - {\partial {u}}/{\partial y}$ and arrows of 
$(u,v)$ components are shown. Line plots show the mode shapes as a function of 
$y$ for 
$\hat {u}$ (——), 
$\hat {v}$ (– – – -) and 
$\hat {w}$ (— - —) components normalized by 
$|\hat {u}|_{max}$.
Figure 6(a–c) show contours of the normalized eigenvalues of the Hessian as a function of 
$(k_x^{+},k_z^{+})$, for increasing observation time 
$t_m^{+}$. Here the Hessian matrices are associated with instantaneous observations of 
$\partial U/\partial y \vert _{wall}$ at 
$t_m^{+}$. In each panel, we identify three important wavenumber pairs, and plot the associated eigenmodes in figure 6(d–f). Modes 
$I_A$ and 
$I_B$ correspond to the two local peaks of the eigenvalues on the horizontal and vertical axis, which represent spanwise and streamwise rolls, respectively. Mode 
$II$ represents large-scale, streamwise-oriented structures with 
$k_z^{+} = 0.02$, or 
$2 {\rm \pi}/ k_z^{+} \approx O(300)$. Compared with other eigenfunctions, mode 
$II$ demonstrates finite sensitivity of wall measurement to the flow beyond the buffer layer towards the core of the channel when the measurement time 
$t_m^{+}$ is large (figure 6f).
Figure 6. The logarithm of the largest eigenvalue of the Hessian matrix for observing 
$\partial U/\partial y \vert _{wall}$ instantaneously at (a–c) 
$t_m^{+} = \{4, 8, 20\}$. Colour and line contours are for 
$Re_{\tau } = 180$ and 
$590$, respectively. The eigenvalues are reported as a function of the horizontal wavenumber vector 
$(k_x^{+}, k_z^{+})$ and are normalized by the respective supremum. (d–f) Profiles of the eigenfunction 
$|\hat {\boldsymbol {u}}|$ associated with the wavenumbers marked in panels (a–c). Thin red lines correspond to 
$Re_{\tau } = 180$ and thick black lines are for 
$Re_{\tau } = 590$.
Five remarks are notable in connection with figure 6. (i) The supremum of the eigenvalues for observing 
$\partial U/\partial y \vert _{wall}$ is associated with a two-dimensional mode 
$I_A$ (
$k_z=0$), which has the form of a roll in the 
$x$–
$y$ plane for all the considered times. (ii) At short time, two perturbations that can effectively influence the wall measurement are modes 
$I_A$ and 
$I_B$. They represent an immediate sensitivity that is not significantly affected by the underlying dynamics of the flow, and their eigenvectors are clustered near the wall. (iii) At longer times, since fine-scale near-wall structures would be dampened by viscous effects, modes 
$I_A$ and 
$I_B$ shift towards longer streamwise and spanwise wavelengths. These larger structures survive longer times, and can be amplified by the flow dynamics over the time duration between initial condition and measurement. (iv) All modes extend farther from the wall at longer observation times, as the domain of dependence of the measurement expands in backward time (cf. figure 2). However, most eigenfunctions vanish above the buffer layer – a detail that we quantify in Appendix A. The practical implication is that wall observations of the streamwise stress are insensitive to many of the turbulence scales above the buffer layer, especially large values of 
$(k_x^{+}, k_z^{+})$. (v) A very important exception is in a small region near mode 
$II$ where the eigenfunctions maintain a finite value in the core of the channel; these modes correspond to the sensitivity of observations of the streamwise wall stress to outer large-scale structures. These elongated motions are immune to the sheltering effect of the strong near-wall shear (see Hunt & Durbin (Reference Hunt and Durbin1999) and Zaki & Saha (Reference Zaki and Saha2009) for discussion of shear sheltering), and their influence on the wall stress is consistent with earlier efforts focused on the forward evolution of turbulence (Abe et al. Reference Abe, Kawamura and Choi2004; Mathis et al. Reference Mathis, Hutchins and Marusic2009; Hwang et al. Reference Hwang, Lee, Sung and Zaki2016; You & Zaki Reference You and Zaki2019). The implication of last two observations (iv) and (v) is evident in figure 1, where only the large-scale motions in the logarithmic layer are accurately reconstructed and the finer scales of turbulence are absent from the estimated state.
The above discussion of the behaviour of the Hessian eigenvalues at long times remains qualitatively unchanged when considering wall observations of instantaneous spanwise shear stress or pressure. The key distinctions are made by aid of figure 7, where the associated largest eigenvalues are plotted as a function of 
$(k_x^{+}, k_z^{+})$, at observation time 
$t_m^{+} = 20$. The most effective flow structure for observing 
$\partial W/\partial y \vert _{wall}$ is three-dimensional with wavenumber 
$(k_x^{+}, k_z^{+}) \approx (0.05, 0.05)$. As for observing the wall pressure, the eigenvalue spectrum shows sensitivity to spanwise rolls and much weaker sensitivity to streamwise ones, the latter being known to develop without an associated strong pressure perturbation (Phillips Reference Phillips1969). The profiles of modes 
$I$ and 
$II$ shown in figure 7 are, respectively, similar to those of modes 
$I_A$ and 
$II$ in figure 6(f) for observing 
$\partial U/\partial y \vert _{wall}$, also at 
$t_m^{+}=20$. Most importantly, for all three types of observations mode 
$II$ maintains a finite value beyond the buffer layer, which represents the sensitivity of wall data to the outer large-scale motions. Mathematically, the eigenmodes of the Hessian matrix can rigorously be interpreted as (a) the leading search directions in adjoint state estimation problem or (b) as the perturbations that yield the largest influence on wall observations. Such interpretations, while accurate, do not guarantee that such eigenstructures are present, or observable, in developed wall turbulence. In addition, we argued on physical grounds that short-time observations may not be sensitive to the dynamics that lead to the generation of energetic flow structures in wall turbulence because such structures develop on longer time scale. We also argued on physical grounds that longer-time observations are sensitive to initial disturbances that amplify by the flow dynamics thus having a large wall signature. In order to support our interpretation, we can tailor the Hessian analysis to focus on the most energetic modes of developed channel flow at different 
$(k_x^{+}, k_z^{+})$. Specifically, we evaluate the sensitivity of wall observations to flow structures 
$\boldsymbol {v}$ obtained from a POD of channel-flow turbulence (Lumley Reference Lumley1967; Moin & Moser Reference Moin and Moser1989; Taira et al. Reference Taira, Brunton, Dawson, Rowley, Colonius, McKeon, Schmidt, Gordeyev, Theofilis and Ukeiley2017).
Figure 7. Contours of the logarithm of the largest eigenvalue of the Hessian matrix for observing (a) 
$\partial W/\partial y \vert _{wall}$ and (b) 
$P\vert _{wall}$, instantaneously at 
$t_m^{+} = 20$, as a function of the horizontal wavenumber vector and normalized by the respective suprema. Colour and line contours correspond to 
$Re_{\tau } = 180$ and 
$590$, respectively. Line plots are profiles of the eigenfunctions 
$|\hat {\boldsymbol {u}}|$ associated with the marked eigenvalues 
$I$ and 
$II$, at 
$Re_{\tau } = 590$.
Assume that the initial perturbation 
$\boldsymbol u_0$ is aligned with a POD mode 
$\boldsymbol u_0 = \alpha \boldsymbol {v}$, where 
$\alpha$ is the amplitude. The cost function (2.12) can be written as
\begin{equation} \mathcal{J}(\alpha \boldsymbol {v};t_m) = \frac{1}{2S}\int_S [\alpha \boldsymbol {v}, \boldsymbol{u}^{*}]^{2}{\rm d}x_m \,{\rm d}z_m = \frac{1}{2S}\int_S \alpha^{2} [ \boldsymbol {v}, \boldsymbol{u}^{*} ]^{2} {\rm d}x_m\,{\rm d}z_m. \end{equation}
Since the POD mode generally satisfies 
$[\boldsymbol {v},\boldsymbol {v}] = 1$, the inner product 
$[\boldsymbol v,\boldsymbol u^{*}]$ is the scalar projection of the adjoint field onto the POD mode. Given 
$\boldsymbol {v}$, the cost function only depends on 
$\alpha$, and the corresponding projected Hessian is the second-order derivative with respect to 
$\alpha$,
\begin{equation} \boldsymbol{\mathcal{H}}_{POD}(t_m) = \frac{\partial^{2} \mathcal{J}}{\partial \alpha \partial \alpha} = \frac{1}{S} \int_S \left[ \boldsymbol {v}, \boldsymbol{u}^{*} \right]^{2} {\rm d} x_m\,{\rm d} z_m. \end{equation}
Note that the projected Hessian is a scalar that quantifies the curvature along a POD direction. This scalar becomes equivalent to the eigenvalues from the preceding section when the POD mode is aligned with a single eigenfunction of the original Hessian. For projection onto POD modes at different 
$(k_x^{+},k_z^{+})$, the corresponding Hessian is denoted as 
$\hat {\boldsymbol {\mathcal {H}}}_{POD}(k_x^{+},k_z^{+};t_m)$.
Analogous to the earlier interpretation, the projected Hessian 
$\hat {\boldsymbol {\mathcal {H}}}_{POD}$ quantify the difficulty of turbulence estimation, but in terms of reconstructing the most energetic structures, or POD modes, in turbulent channel flow. Larger values imply a significant alignment between the eigenmodes of the original Hessian and the POD modes, and the associated flow structures will therefore be effectively reconstructed by adjoint-variational estimation.
Results are reported in figure 8 for 
$Re_{\tau } = 180$. The colour contours are the eigenvalues 
$\lambda ^{\prime }$ of the Hessian as a function of 
$(k_x^{+},k_z^{+})$, when the projection is onto (figure 8a, left subpanels) the first and (figure 8a, right subpanels) the second POD modes; for comparison, the eigenvalues of the original Hessian, without projection, are included as lines. The maximum eigenvalues of the latter are used for normalization. For early observations, 
$t_m^{+}=4$, the spectra of the original and projected Hessian are dissimilar: specifically, the original peak eigenvalue at 
$k_x^{+} = 0.12$ (mode 
$I_A$) does not persist, which indicates that the wall observations are not sensitive to the energetic turbulent structures in that wavenumber range. Instead, the sensitivity at 
$k_x^{+} = 0.12$ (mode 
$I_A$) is to the near-wall, small scales; the adjoint-variational data assimilation procedure is therefore effective at reconstructing these structures. For late observations (
$t_m^{+}=20$), the contours of the spectra of both the original and projected Hessian are more similar near the peak 
$(k_x^{+}, k_z^{+}) = (0, 0.08)$, which is marked ‘A’ in figure 8 and is associated with streamwise-elongated structures with spanwise width 
$2 {\rm \pi}/ k_z^{+} \sim O(100)$. The implication is that the maximum sensitivity of wall observations is aligned with the most energetic POD modes at this wavenumber. Adjoint-variational data assimilation using observations at 
$t_m^{+}=20$ will therefore target a reconstruction of these energetic structures in the initial state, because the underlying flow dynamics amplify these structures in forward time and thus they lead to the largest impact on the observations.
Figure 8. (a) Eigenvalues of the projected Hessian showing the sensitivity of wall observations of 
${\partial U}/{\partial y}$ to (left subpanels) the first and (right subpanels) the second POD modes of turbulent channel flow at 
$Re_{\tau } = 180$. Lines are reproduced from figure 6, for the original Hessian without projection. Observation times are (top subpanels) 
$t_m^{+}=4$ and (bottom subpanels) 
$t_m^{+}=20$. (b) The first 
$\hat {\boldsymbol {v}}_1$ and second 
$\hat {\boldsymbol {v}}_2$ POD modes for wavenumbers (‘A’, solid) 
$(k_x^{+}, k_z^{+}) = (0, 0.08)$ and (‘B’, dashed) 
$(k_x^{+}, k_z^{+}) = (0, 0.02)$. (c) Eigenmodes of the original Hessian 
$|\hat {\boldsymbol {u}}^{*}|$ for the same wavenumber pairs, at 
$t_m^{+} = 20$.
A more detailed interpretation of the spectra is aided by comparison of the leading POD modes (figure 8b) and the eigenfunctions of the original Hessian (figure 8c). The results for point ‘A’ show that the peak projection at 
$(k_x^{+}, k_z^{+}) = (0, 0.08)$ is due to the similarity between the modes below the buffer layer, both being associated with streamwise-elongated and energetic near-wall streaks. Another wavenumber pair of interest is marked ‘B’ in figure 8, at 
$(k_x^{+}, k_z^{+}) = (0, 0.02)$ which was previously labelled mode 
$II$ in the analysis of the original Hessian. The eigenvalues at ‘B’ become appreciable when 
$t_m^{+}=20$. The POD modes (figure 8b) are finite throughout the majority of the channel height, and the eigenfunction of the original Hessian (figure 8c, also mode 
$II$ in figure 6f) is similarly finite outside the buffer layer at this 
$t_m^{+}$. The projected Hessian at ‘B’ thus satisfies two important conditions: (i) the eigenvalue 
$\lambda ^{\prime }_1$ is appreciable and (ii) both the POD and original-Hessian modes are aligned in the outer flow. The practical implication is that wall observations are sensitive to the outer streamwise-aligned (
$k_x^{+}=0$) large-scale (
$2{\rm \pi} /k_z^{+} \sim O(300)$) energetic structures, and 4DVar can estimate these structures in the initial state from wall measurements.
The above discussion was framed in terms of the eigenspectrum of the Hessian, and the reported eigenvalues were normalized in all the figures by their supremum at the corresponding observation time. The time dependence of the supremum is reported in Appendix A, and shows long-time exponential amplification due to the Lyapunov behaviour of the adjoint system. As a result, when observations are accumulated during a long time horizon, for example 
$t_m^{+} = 50$ as in § 2.1, the estimation of the initial state is most affected by the late observations. It is instructive to compare the notions of Lyapunov divergence in the adjoint and forward Navier–Stokes equations: just as in the forward problem where an initial perturbation to the state grows exponentially in time with the Lyapunov exponent, an initial perturbation in the adjoint equations amplifies exponentially in backward time. For example, at 
$Re_\tau = 180$ Nikitin (Reference Nikitin2018) reported a Lyapunov exponent 
$\sigma ^{+}=0.021$ for the forward problem; the adjoint field has the same exponential rate as we will report in § 3.3. The interpretation of the adjoint behaviour is, however, different from the familiar forward problem. In the context of the Hessian at optimality, the sensitivity of an observation, which is represented by the adjoint field starting from the observation kernel, leads to exponentially ‘larger Hessian’ at 
$\tau = t_m$. Not only does the supremum eigenvalue increase, but so does its separation from smaller ones as well (cf. figure 6). As a result, the condition number of the Hessian becomes larger and it becomes much more difficult to solve the state-estimation problem, accurately. The important implication is that errors in observations, in particular at larger 
$t_m$, can strongly obscure the reconstruction of the initial state.
3.3. Asymptotic behaviour of the adjoint field and of the gradient of the cost function
The analysis in the previous section focused on evaluation of the Hessian when the gradient of the cost function vanishes. The results highlighted the difficulty of accurately predicting the full turbulent state, especially beyond the buffer layer, from wall observations. Even when the initial estimate of the flow is near the true state, the wall stress has a diminishing sensitivity to the state with distance from the wall and, in the outer region, can improve the estimation of the large-scale motions only. Generally, however, an initial guess of the state may just be the mean-flow profile that is far from the true field. In addition, the assimilation window of observations is often long. As such, the exponential amplification of the adjoint in backward time is evidence of a very large gradient of the cost function with respect to the initial condition, as well as an ill-conditioned Hessian. These realities impose severe restrictions on the step size in the gradient descent method. In this section, we examine the kinetic energy of the adjoint field because it is directly proportional to the magnitude of the gradient of the cost function with respect to the initial condition. The results will highlight regions of the flow that contribute most to this amplification, and will be contrasted to the evolution of perturbations in the forward field for which we have established understanding.
In § 3.1, we examined the adjoint fields due to impulses at the wall (figures 2 and 3), at short observation times. At long reverse times, the initially localized patches spread, become chaotic and fill the domain. Ultimately they reach a statistical state that is independent of the observation type and location. Therefore, the long-time asymptotic behaviours of the adjoint fields from (2.5) and (2.11) are the same; the former provides the gradient of the cost function and the latter is associated with Hessian matrix. The asymptotic statistical state can be reached efficiently by performing adjoint computations starting from broadband adjoint velocities throughout the computational domain, which is the approach adopted here. It should be noted that the adjoint equations are linear, and hence the initial field is anticipated to grow exponentially in backward time, indefinitely. This behaviour can be likened to the amplification of perturbations to a forward system: while the nonlinear evolution of the perturbations leads to a statistically saturated state, when the linearized Navier–Stokes equations are adopted the perturbations amplify exponentially, unabated.
Starting from initial random perturbations that were projected onto solenoidal fields, we simulated the flow evolution using both the linearized forward (2.7) and adjoint (2.11) Navier–Stokes equations in order to contrast the long-time statistical properties of both models. Recall that in both cases the base state in the governing equations is the three-dimensional, time-dependent, fully turbulent flow at the corresponding Reynolds number. Figure 9(a) shows the evolution of the perturbation energy,
\begin{equation} K(t) = \frac{1}{2V}\int_V | \boldsymbol{u} |^{2} \,{\rm d}V \quad\text{and}\quad K^{*}(t) = \frac{1}{2V}\int_V | \boldsymbol{u}^{*} |^{2} \,{\rm d}V. \end{equation}
The integration is performed over the full three-dimensional domain, and the exponential growth rate 
$\sigma ^{+} \approx 2.11 \times 10^{-2}$ at 
$Re_{\tau } = 180$ agrees with the value reported by Nikitin (Reference Nikitin2018) for forward evolution.
Figure 9. (a) Exponential energy amplification for the (solid) linearized forward and (dashed) adjoint fields. Increasing line thickness and vertical shift correspond to 
$Re_{\tau } = \{180, 590, 1000\}$. For each Reynolds number, the Lyapunov exponents are similar for the forward (
$\sigma ^{+} = \{2.11, 2.48, 2.73\}\times 10^{-2}$) and adjoint (
$\sigma ^{+} = \{2.04, 2.33, 2.61\}\times 10^{-2}$). (b) Side views of sample (top subpanel) linearized forward and (bottom subpanel) adjoint fields during the statistically stationary state for 
$Re_{\tau } = 180$.
Beyond a short-lived transient, the statistical behaviour of the perturbation field within the exponential regime is unchanged, to within a rescaling. We therefore normalize quantities by their Lyapunov amplification,
\begin{equation} \mathring{\boldsymbol{u}} = \boldsymbol{u} {\rm e}^{-\sigma t}, \quad \mathring{\boldsymbol{u}}^{*} = \boldsymbol{u}^{*} {\rm e}^{-\sigma \tau}, \end{equation}
and analyse 
$\mathring {\boldsymbol {u}}$ and 
$\mathring {\boldsymbol {u}}^{*}$ as statistically stationary fields. Figure 9(b) shows instantaneous side views of the normalized forward and adjoint perturbation fields for 
$Re_{\tau } = 180$. The adjoint field reflects the more localized near-wall behaviour, which is manifest in the form of concentrated high-intensity patches of turbulence. The horizontally and time-averaged kinetic energy of the normalized variables,
\begin{equation} \overline{\mathring{k}^{(*)}}= \tfrac{1}{2}\,\overline{| \mathring{\boldsymbol{u}}^{(*)} |^{2}}, \end{equation}
was evaluated and the wall-normal profiles are plotted in figure 10. Each profile is normalized to unit integral, and is plotted multiplied by 
$y^{+}$ to ensure that the area under the curve is representative of the integral. The energy of the forward perturbations peaks in the buffer layer, and decays within the logarithmic region towards the channel centre. The energy of the adjoint field exhibits narrower support, with a more concentrated peak and much faster decay within the buffer layer. The large gradient of the cost function in the state estimation problem is therefore most concentrated in the buffer layer. These results are also relevant to the poor condition number of the Hessian. In the previous section (§ 3.2), figure 6(a–c) shows that the separation between the supremum and smaller eigenvalues of the Hessian increases with 
$t_m$, and hence the Hessian becomes progressively more ill-conditioned. The asymptotic behaviour of the adjoint system (figures 9–10) reinforces this statement by showing that observation errors could lead to exponential amplification of errors in the buffer layer, due to adjoint chaos. Taken altogether, the large gradient in the buffer layer and ill-conditioned Hessian render the state estimation problem very challenging.
Figure 10. Profiles of the horizontally and time-averaged kinetic energy in the (black solid) linearized forward and (red dashed) adjoint velocity fields. Increasing line thickness corresponds to Reynolds numbers 
$Re_{\tau } = \{180, 590, 1000\}$. Each curve is normalized to have unit integral.
The evolution equations of the kinetic energies of the perturbations are derived by performing the dot product of the linearized forward Navier–Stokes equations (2.7) with 
$\boldsymbol {u}$ and similarly the dot product of the adjoint Navier–Stokes equations (2.11) with 
$\boldsymbol {u}^{*}$. The resulting equations are
\begin{equation} \left.\begin{gathered} \frac{\partial k}{\partial t} ={-}u_i u_j S_{ij} - \frac{\partial }{\partial x_j} \left(\frac{1}{2}u_i u_i U_j + u_j p -\frac{2}{Re} u_i s_{ij}\right) - \frac{2}{Re} s_{ij} s_{ij},\\ \frac{\partial k^{*}}{\partial \tau} ={-}u^{*}_i u^{*}_j S_{ij} - \frac{\partial }{\partial x_j} \left(- \frac{1}{2}u^{*}_i u^{*}_i U_j - u^{*}_j p^{*} -\frac{2}{Re} u^{*}_i s^{*}_{ij}\right) - \frac{2}{Re} s^{*}_{ij} s^{*}_{ij}, \end{gathered}\right\} \end{equation}
where 
$k = (u_i u_i)/{2}$ is the kinetic energy, 
$s_{ij} = (\partial u_i/\partial x_j + \partial u_j / \partial x_i)/{2}$ and 
$S_{ij} = (\partial U_i / \partial x_j + \partial U_j/\partial x_i)/{2}$ are the rate of strain tensors for the perturbation and the reference fields, respectively, and stars denote adjoint quantities. We multiply the forward equation by 
$\exp (-2 \sigma t )$ and the adjoint by 
$\exp (-2 \sigma \tau )$, and average each of them in time and in the homogeneous 
$x$ and 
$z$ directions, which yields
\begin{gather} 2\sigma \bar{\mathring{k}} = \underbrace{- \overline{\mathring{u}_i \mathring{u}_j S_{ij}}}_{\mathcal{P}} \underbrace{- \dfrac{{\rm d}}{{\rm d} y}\left(\dfrac 12 \overline{\mathring{u}_i \mathring{u}_i U_2} +\overline{\mathring{u}_2 \mathring{p}} - \dfrac{2}{Re} \overline{\mathring{u}_i \mathring{s}_{i2}}\right)}_{\mathcal{T}} - \underbrace{\dfrac{2}{Re} \overline{\mathring{s}_{ij}\mathring{s}_{ij}}}_{\epsilon}, \end{gather}
\begin{gather}2\sigma \overline{\mathring{k}^{*}} = \underbrace{- \overline{\mathring{u}^{*}_i \mathring{u}^{*}_j S_{ij}}}_{\mathcal{P}^{*}} \underbrace{- \dfrac{{\rm d}}{{\rm d} y}\left(-\dfrac{1}{2}\overline{\mathring{u}^{*}_i \mathring{u}^{*}_i U_2} - \overline{\mathring{u}^{*}_2 \mathring{p}^{*}} - \dfrac{2}{Re} \overline{\mathring{u}^{*}_i \mathring{s}^{*}_{i2}}\right)}_{ \mathcal{T}^{*}} - \underbrace{\dfrac{2}{Re} \overline{\mathring{s}^{*}_{ij}\mathring{s}^{*}_{ij}}}_{\epsilon^{*}}. \end{gather} The various terms in the kinetic energy equations (3.7) and (3.8) are reported in figure 11. For the forward flow, the profiles are qualitatively similar to the recent results by Nikitin (Reference Nikitin2018) at 
$Re_{\tau }$ around 
$390$. In contrast to the forward terms, the adjoint profiles are confined to a narrower wall-normal extent 
$y^{+} \le 60$, have more compact support and have a relatively more prominent peak in the buffer layer where the adjoint activity is most intense. The cause of the difference between the forward and adjoint terms can perhaps be understood in terms of the reversal in the turbulent motions in the adjoint evolution: sweeps that are prominent near the wall in the forward dynamics become ejections in the adjoint, and pronounced ejections above the buffer layer in the forward equations become sweeps in the adjoint. Ultimately the very large concentrated production of the adjoint energy in the buffer layer is responsible for the exponential amplification in reverse time. In terms of the state-estimation problem, this corresponds to extreme gradients in the cost function and an ill-conditioned Hessian with expanding assimilation horizon; the buffer layer thus obscures the interpretability of wall observations. However, when these results are viewed through the lens of the spectral analysis of the Hessian (§ 3.2), it becomes evident that wall observations remain sensitive to low 
$k_x$ flow structures beyond the buffer layer, which have been related to outer–inner interactions in wall turbulence. Taken all together, the present statistical results and the Hessian analysis provide a comprehensive understanding of the interpretability of wall observations in the attempt to estimate the initial condition of channel flow turbulence.
Figure 11. Kinetic energy budget for (a) the linearized forward and (b) the adjoint equations. Results for 
$Re_{\tau } = \{180, 590, 1000\}$ and shown with increasing line thickness. Curves are normalized such that 
$\max _y 2\sigma \bar {\mathring {k}}=\max _y 2\sigma \overline {\mathring {k}^{*}} = 1$.
4. Conclusion
The nonlinear and chaotic natures of turbulence render its reconstruction from limited observations a challenging problem. The present work focused on the canonical configuration of turbulent channel flow. While all the vorticity in this flow has its origin at the wall, estimating the turbulent state from observations of the wall shear stresses of wall pressure is notoriously difficult.
The present effort focuses on the interpretability of the wall measurements, from the perspective of adjoint-variational data assimilation (i.e. 4DVar). In this approach, the estimation problem is formulated as a constrained optimization where we attempt to identify the initial condition whose nonlinear Navier–Stokes evolution reproduces all the available observations. Discrepancies between the model predictions and the available observations define the cost function to minimize, and which also features in the forcing term to the adjoint equations that are marched back in time. The outcome of one forward-adjoint loop is the gradient of the cost function that is used in the gradient-based minimization procedure. The accuracy of the state estimation depends on a number of factors, including the sensitivity of the observations to the flow state which is related to the geometry of the cost function, most importantly its gradient and Hessian. In order to ensure accuracy of our computations, a discrete adjoint is adopted which satisfies the forward-adjoint duality relation to machine precision.
In order to frame the discussion, we provided a summary of recent results by Wang & Zaki (Reference Wang and Zaki2021) for estimating the initial state of turbulent channel flow from wall observations, at friction Reynolds numbers 
$Re_\tau = \{100, 180, 392, 590\}$. The assimilated states were obtained starting with an initial guess from an LSE and performing 100 iterations of the 4DVar algorithm, for each of the reported Reynolds numbers. The predicted flow state displayed important characteristics. (i) The near-wall turbulence was accurately reconstructed with a nearly perfect correlation with the true flow that generated the observations. This region therefore diminishes in physical height at higher Reynolds numbers. (ii) The correlation coefficient between the estimated and true states, however, decays precipitously across the buffer layer and into the channel core. Despite the reduction in accuracy, the outer large-scale energetic motions which modulate the near-wall flow are captured by the estimation.
The above characteristics were explained in terms of the domain of dependence of observations, and by examining geometric properties of the cost function for state estimation from instantaneous wall measurements. Specifically, we considered a scenario where the estimate of the initial flow field is infinitesimally close to the true solution and an instantaneous observation at 
$t=t_m$ is used to improve the estimate. At optimality the gradient of the cost function vanishes and its Hessian determines the behaviour of the optimization procedure. In order to evaluate the Hessian, an efficient approach is introduced that exploits a forward-adjoint duality relation. The relevant adjoint field in this relation is initiated from an observation kernel at the wall; its spatiotemporal evolution in reverse time is the domain of dependence of the kernel, which has a universal behaviour across Reynolds numbers when scaled in viscous units. The interpretation of this adjoint field as the domain of dependence of observations can be exploited in the design of measurements, for example to optimally distribute the measurements and ensure sensitivity to the entire state. The Hessian matrix is the cross-correlation of this adjoint field at the final reverse time, 
$\tau = t_m$. It was evaluated efficiently in Fourier space for each horizontal wavenumber pair 
$(k_x,k_z)$, and was integrated over observation locations 
$(x_m,z_m)$ in the plane of the wall. The eigenmodes of the Hessian were then analysed in order to demonstrate the capacity to use wall observations for 4DVar estimation of the state of channel-flow turbulence.
When observing the streamwise wall shear-stress at time 
$t^{+}_m = 4$ (figure 6), the largest Hessian eigenvalue, or curvature, is for mode 
$I_A$ with 
$(k_x^{+}, k_z^{+}) \approx (0.12, 0)$ and the associated eigenvector is a spanwise roll. The measurements are therefore most sensitive to a disturbance in the initial state that is aligned with mode 
$I_A$ and, from the perspective of state estimation, this component of the initial condition is reliably targeted by the 4DVar procedure. For longer measurement times, similar spanwise rolls with smaller 
$k_x^{+}$ have the greatest influence on the measurement. For most wavenumber pairs, the leading eigenvectors of the Hessian are concentrated near the wall (
$y^{+} \le 40$), and hence the turbulence is accurately estimated in this region. One important exception is mode 
$II$ at 
$(k_x^{+}, k_z^{+}) \approx (0, 0.02)$, which by 
$t_m^{+}=20$ becomes finite outside the buffer layer and represents the sensitivity of wall measurements to outer large-scale structures.
The role of the observation time 
$t_m$ is noteworthy: at small values of 
$t_m^{+} \approx 4$, the Hessian eigenspectrum captures the sensitivity of wall observations to perturbations to the flow state that may not have sufficient time to amplify or decay due to the dynamics of the forward operator. In contrast, at larger values of 
$t_m^{+} \approx 20$, the Hessian eigenspectrum shows higher sensitivity to lower wavenumbers because the associated perturbations to the flow can amplify in forward time and have a sizable impact on wall signature, while high-wavenumbers are dissipated by viscosity. This perspective was reinforced by projecting the Hessian onto the first and second modes from POD of channel flow turbulence (§ 3.2). After projection, the Hessian yields the curvature of the cost function along the directions of these POD structures in channel flow, or the capacity of state estimation to predict them. Another equally valid interpretation of the projected Hessian is the sensitivity of wall observations to the energetic POD structures. The peak sensitivity of the original Hessian at (mode 
$I_A$) is not retained in the projection because it was not associated with energetic flow structures but rather with near-wall small-scale turbulence. At larger times 
$t^{+}_m = 20$, the contours of the eigenvalues of the original and projected Hessians are more aligned, both showing strong sensitivity of wall observations to streamwise-elongated and energetic flow structures.
The increase in observation time is accompanied by another important effect: the largest eigenvalues of the Hessian increase exponentially and their separation from smaller ones becomes larger (cf. figure 6), which renders the Hessian progressively more ill-conditioned. If observations are accumulated during a time horizon, late observations dominate the state estimation problem and the solution of the inverse problem becomes prone to errors. Specifically, infinitesimal uncertainties in later observations can lead to large errors in the estimation of the initial state. The condition number of the Hessian also provides an objective approach to optimize the design of measurements, for example by optimally placing sensors or weighting their relative contributions to the cost function in order to minimize the condition number.
The starting estimate for practical data assimilation often deviates appreciably from the true flow state, and hence the gradient of the cost function as well as its Hessian are important. Both quantities are related to the adjoint field, and are affected by adjoint chaos at long reverse times which is markedly different from the more familiar forward counterpart. We contrasted the evolution of small-amplitude perturbations to the forward flow field and the adjoint model, using the linearized Navier–Stokes equations and their adjoint, respectively. The long-term statistical behaviour of the adjoint chaos in turbulent channel flow was reported for the first time. The results provided a glimpse of the cause for the exponential growth of the adjoint in backward time. The energy of the adjoint variable peaks in the buffer layer, and its distribution has narrower support in the wall-normal direction than that of the forward variable. By rescaling the energy equation by twice the Lyapunov amplification rate of perturbations, we are able to compute the different contributions to its budget. Production is concentrated in the buffer layer and decays quickly away from its peak, relative to the profile for the forward problem. This high adjoint turbulent kinetic energy of its production is the root of chaos in the buffer layer. These adjoint characteristics are indicative of the extreme gradients of the cost function and its ill-conditioned Hessian matrix, which both obscure the reconstruction of the initial state from wall data. These results support the practice of limiting the assimilation time horizon based on the Lyapunov time scale (Li et al. Reference Li, Zhang, Dong and Abdullah2020; Wang & Zaki Reference Wang and Zaki2021), in order to mitigate the impact of adjoint (and also forward) chaos, and adopting a sliding window (Chandramouli, Mémin & Heitz Reference Chandramouli, Mémin and Heitz2020) or a cycling scheme (Fisher & Auvinen Reference Fisher and Auvinen2011) for longer horizons. In addition, while we can precisely reconstruct all the turbulence scales in the near-wall region from boundary observations, the buffer layer acts as a low-pass filter beyond which we were able to accurately estimate only the low-frequency, large-scale, energetic outer motions.
Acknowledgements
Computational resources were provided by the Maryland Advanced Research Computing Center (MARCC).
Funding
The authors acknowledge financial support from the Office of Naval Research (grant N00014-20-1-2715).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Temporal behaviour of the Hessian eigenmodes
The eigenmodes of the Hessian represent the sensitivity of the observations to perturbations to the flow state. The results in the main text showed that, for the majority of modes, the leading eigenfunctions are concentrated in the near-wall region, expand vertically with observation time and decay above the buffer layer (see figure 6). One important exception was noted for the modes with horizontal wavenumber vector 
$(k_x^{+}, k_z^{+}) \approx (0, 0.02)$ at 
$t_m^{+} \ge 20$ that provide a connection between wall observations and the outer large-scale structures. In order to characterize the general wall-normal extent of penetration of the Hessian eigenmodes in the channel, we adopt two measures, 
$y_c$ and 
$y_{0.1}$. The first is the moment
\begin{equation} y_c \equiv \frac{\displaystyle\int_0^{1} y |\hat{\boldsymbol{u}}| \,{{\rm d} y}}{ \displaystyle\int_0^{1} |\hat{\boldsymbol{u}}| \,{{\rm d} y}}, \end{equation}
where 
$|\hat {\boldsymbol {u}}| \equiv \sqrt {|\hat {u}|^{2}+|\hat {w}|^{2}+|\hat {w}|^{2}}$, and the second is the height at which the eigenvector decays to 10 % of its maximum value,
\begin{equation} y_{0.1} = \max{ \left\{0< y<1 ~|~\left({|\hat{\boldsymbol{u}}(y)|}/{|\hat{\boldsymbol{u}}|_{max}}\right) = 0.1\right\}}. \end{equation} For the Reynolds number of 
$Re_{\tau } = 180$, we considered the Hessian eigenspectrum due to an instantaneous observation at time 
$t_m$, and averaged 
$y_c$ and 
$y_{0.1}$ over all horizontal wavenumber vectors, 
$\overline {y_c} \equiv \langle y_c\rangle _{k_xk_z}$ and 
$\overline {y_{0.1}} \equiv \langle y_{0.1}\rangle _{k_xk_z}$. The results are shown as a function of 
$t_m^{+}$ in figure 12. While both measures initially increase with observation time, they level off beyond a few time units. The value of 
$\overline {y_c}$ never exceeds the buffer layer, while the less conservative estimate remains in the range 
$\overline {y_{0.1}^{+}} < 50$. The elongated modes which penetrate beyond the buffer layer are therefore not reflected in the average, but are important to recall since they enable the reconstruction of the outer large-scale structures.
Figure 12. Measures of persistence of the Hessian eigenfunctions away from the wall, 
$\bar {y}_c$ and 
$\bar {y}_{0.1}$. These quantities are evaluated for the leading eigenvector at every horizontal wavenumber pair, and then averaged over all 
$(k_x^{+}, k_z^{+})$. The Reynolds number is 
$Re_{\tau } = 180$. Results correspond to observing (——) 
$\partial U/\partial y \vert _{wall}$, (– – – -) 
$\partial W/\partial y \vert _{wall}$ and (— - —) 
$P \vert _{wall}$. (a) Observation at one time instance 
$t_m$. (b) Observations at every time instance within the assimilation window 
$0 < t^{\prime } \le t_m$.
Since in practice observations are often available not just at one time instance but rather over an assimilation time window, the cost function is integrated in time, i.e. all the available observations are used to estimate of the initial condition. The relevant cost function is therefore the time integrated form. For example, when we observe 
$\partial U/\partial y \vert _{wall}$, the time-integrated cost function is
\begin{equation} J_{u} = \frac{1}{2S}\int_0^{t_m}\int_S \left( \left.\frac{\partial u}{\partial y}\right\vert_{(\boldsymbol{x}_m,t^{\prime})} \right)^{2} \,{\rm d}x_m \,{\rm d}z_m \,{\rm d}t^{\prime} \end{equation}and the associated Hessian involves the time integral
\begin{equation} H_u = \frac{\partial^{2} J_u}{\partial \boldsymbol{u}_0\partial \boldsymbol{u}_0} =\frac{1}{S} \int_0^{t_m}\int_S \boldsymbol{u}^{*}\left({\bullet};\boldsymbol{x}_m,t^{\prime}\right)\boldsymbol{u}^{*}\left({\bullet};\boldsymbol{x}_m,t^{\prime}\right) {\rm d}x_m \,{\rm d}z_m \,{\rm d}t^{\prime}, \end{equation}
where 
$\bullet$ refers to 
$(\boldsymbol x, t=0)$. We can similarly evaluate 
$H_w$ and 
$H_p$. In other words, the Hessian associated with estimating the initial state from the full time-history of observations is the time integral of the Hessian associated with isolated, or instantaneous, observations. The integral is dominated by late observations because the adjoint field amplifies exponentially with 
$t_m$. As such, we can anticipate that the results from figure 12(a) remain qualitatively unchanged when observations are accumulated over the assimilation window. This expectation is verified in figure 12(b) where 
$\overline {y_c}$ and 
$\overline {y_{0.1}}$ are plotted, evaluated from the eigenfunction of the time-integrated Hessian. These results supplement the discussion of figure 6 in the main text, and provide statistical confirmation of the challenge of estimating channel-flow turbulence from wall observations.
The discussion at the end of § 3.2 referred to the time-dependence of the supremum eigenvalue of the Hessian associated with instantaneous observations. This dependence is reported in figure 13 as a function of 
$t_m^{+}$. After an initial transient, the eigenvalues increase exponentially due to the chaotic nature of the adjoint system. This trend, and the larger separation between the supremum and infimum (figure 6) with observation time, demonstrate that the Hessian becomes increasingly ill-conditioned. As a result, the state-estimation problem becomes progressively more challenging.
Figure 13. The supremum eigenvalue of the Hessian matrices associated with instantaneous observations at all horizontal wavenumber pairs 
$(k_x^{+}, k_z^{+})$ as a function of the observation time 
$t_m^{+}$ for 
$Re_{\tau } = 180$. Different wall-observation modalities are shown: 
$\partial U/\partial y\rvert _{wall}$ (——); 
$\partial W/\partial y\rvert _{wall}$ (– – – -); and 
$P\rvert _{wall}$ (— - —). Lyapunov amplification at 
$\sigma ^{+} = 0.0215$ is marked for comparison.
Open access
