State estimation in turbulent channel flow from limited observations

The flow evolution is governed by the incompressible Navier-Stokes equations,

where $t$ is time and $p$ is pressure. These equations are also referred to as the forward model because they are adopted to advance the flow state in time.

The Navier-Stokes equations are solved using a fractional-step method with a local volume flux formulation on a staggered grid (Rosenfeld et al., 1991). The advection terms are discretized by the Adams-Bashforth scheme, and the Crank-Nicolson scheme is adopted for the diffusion terms. The pressure Poisson equation is solved using Fourier transform in the periodic directions and tri-diagonal inversion in wall-normal direction. The algorithm has been applied in a number of direct numerical simulations of transitional and turbulent flows (Zaki et al., 2010; Zaki, 2013; Zaki & Durbin, 2005; Lee & Zaki, 2017). For simplicity, the discretized Navier-Stokes equations will be denoted as,

where $\mathbf{q}^{n}$ is the state vector, including the velocity and pressure at every grid points, and $\mathbf{G}^{n}$ is a matrix that represents the discretized Navier-Stokes operator. Note that $\mathbf{G}^{n}$ is also a function of $\mathbf{q}^{n}$ because the equations are nonlinear.

The true state is statistically stationary turbulence, and is sustained by a known constant pressure gradient in the streamwise direction. Except in §3.5 where we explore the effect of Reynolds numbers, we set $Re_{\tau}=180$. While the forward model at these conditions has been extensively studied (Kim et al., 1987; Jelly et al., 2014), this Reynolds number is higher than previously attempted in the context of adjoint-variational state estimation in channel flow. The domain size and the grid resolution are summarized in table 1. The computational domain is the same as the one adopted by Kim et al. (1987) who used a pseudo-spectral algorithm. For our finite-volume scheme, we have doubled the resolution in each direction and performed extensive validation (see e.g. Jelly et al., 2014). The grid resolution is also reported in viscous units, denoted by superscript $(\cdot)^{+}$, $\Delta x^{+}\equiv(\Delta x/h)Re_{\tau}$. The time step size is $\Delta t^{+}\equiv(\Delta tU_{b}/h)(Re_{\tau}^{2}/Re)=0.058$ such that the Courant–Friedrichs–Lewy (CFL) number is lower than half.

We consider two types of observations: sub-sampled velocity data and stresses on both channel walls. The observations setup is summarized in table 2. In all cases, the estimation window is $T=4.5$ ($T^{+}=50$). This choice is motivated by the following considerations: The duration $T$ should be sufficiently long that each point in the fluid is within the domain of dependence of observations. It should also be longer than the time to “fill” the turbulence energy spectra starting solely form observations. Finally, $T$ should not appreciably exceed the Lyapunov timescale ($\tau^{+}_{\sigma}=48$ at $Re_{\tau}=180$ according to Nikitin, 2018). If $T\gg\tau_{\sigma}$, any infinitesimal perturbation in the initial condition will exponentially amplify and the accuracy of the state estimation deteriorates (Li et al., 2020; Chandramouli et al., 2020). We start with a benchmark case (§3.1, case B1), where the velocity field is observed every eighth point in each dimension, including space and time. The velocity at the observation locations is assumed to be known precisely, without measurement noise. Subsequently (§3.2, cases N1-N2), the data will be contaminated with Gaussian random noise with standard deviation proportional to the local velocity, for example for the streamwise velocity,

The effect of spatial and temporal data resolution is explored in §3.3 (cases RZ1$-$RZ4, RX1$-$RX4, RT1$-$RT4), and in §3.4 the velocity observations in the viscous sublayer and buffer layer are removed. Finally, in §3.5, we evaluate the influence of Reynolds number when the only observations are the wall stresses.

In order to minimize the cost function (1) while satisfying the Navier-Stokes constraint, we introduce the Lagrangian,

Note that the Lagrangian is a function of $\left\{\mathbf{q}\right\}_{n=0}^{N}$ and $\left\{\mathbf{q}^{{\dagger}}\right\}_{n=1}^{N}$. Taking the derivative of the Lagrangian with respect to $\mathbf{q}^{\dagger(n)}$ and setting it to zero, we obtain the forward Navier-Stokes equations (4). By setting the derivatives of the Lagrangian with respect to $\mathbf{q}^{n}$ to zero ($1\leq n\leq N-1$), we obtain the discrete adjoint equations,

which are marched backward in time, and are forced by the gradient of the cost function with respect to the state. The forward operator $\mathbf{G}^{n}$ contains the forward variables $\mathbf{q}^{n}$, which means that the full spatiotemporal evolution of the forward fields are required and must be stored for the solution of the adjoint equations; The second term on the right-hand side of (7) can be analytically derived from the expression of the cost function. When the adjoint equations are marched back to $n=0$, the following relation is obtained,

The initial adjoint field is therefore the gradient of the Lagrangian, and also the gradient of the cost function when both the forward and adjoint equations are satisfied,

Note that $\mathbf{q}^{0}$ contains the velocity field $\mathbf{u}^{0}$ only, because the initial pressure field is not required to solve the incompressible Navier-Stokes equations. Similarly, the initial adjoint field $\mathbf{q}^{{\dagger}0}$ is comprised of $\mathbf{u}^{{\dagger}0}$only. Since the above derivation starts from the discrete Navier-Stokes equations, the gradient obtained using the discrete adjoint in equation (8) is accurate to machine precision. Detailed expressions of the discrete adjoint and verification the forward-adjoint relation are provided in (Wang et al., 2019a).

With the gradient of the cost function, we adopt the Limited-memory Broyden–Fletche–Goldfar–Shanno (L-BFGS) optimization algorithm to minimize the cost function (Nocedal, 1980). In order to guarantee that the estimated initial condition is divergence free, we slightly modify the L-BFGS algorithm by introducing a symmetric projector. The basic idea is to update the new estimate of the initial condition using,

where subscript $k$ denotes the $k$-th iteration of the optimization procedure, and $\mathbf{d}_{k}=-\mathbf{B}_{k}\mathbf{P}\mathbf{u}^{{\dagger}0}$ is the search direction; the matrix $\mathbf{B}_{k}$ is a rank-2 approximation of the inversion of the Hessian matrix of the cost function; and the matrix $\mathbf{P}$ is a symmetric projector, which projects any velocity field $\mathbf{u}^{0}$ or gradient $\mathbf{u}^{{\dagger}0}$ onto the divergence-free space. The symmetry of $\mathbf{P}$ ensures that $\mathbf{P}\mathbf{u}^{{\dagger}0}$ is the gradient of the cost function with respect to $\mathbf{u}^{0}$ when $\mathbf{u}^{0}$ is projected and before being advanced by the forward equations. The step size $\alpha_{k}$ is computed by the line search routine CVRSCH (Moré & Thuente, 1994), which enforces the strong Wolfe condition and adopts cubic interpolation to update $\alpha_{k}$.

Combining the forward solver (§2.1) with the adjoint solver and optimization algorithm (§2.2), we obtain the adjoint-variational state estimation algorithm. A summary is provided in Algorithm 1. In all the examined configurations, the algorithm is always performed for one hundred L-BFGS iterations, and as such comparisons are made using the same computational cost. It is important to emphasize that, in addition to the storage requirements associated with saving the forward fields at full spatiotemporal resolution, the computational cost is also substantial because each of the one hundred L-BFGS iterations comprises at least one forward and one adjoint computations.

In order to provide a qualitative perspective on the performance of algorithm, figure 3 shows realizations of the flow from case B1 (see table 2), at both the (i) initial and (ii) final ($t^{+}=T^{+}=50$) times within the assimilation window. At each instant, the field is visualized using ($a$) the observations data only, ($b$) the adjoint-variational prediction, and ($c$) a detailed comparison of the predicted (colour) and true (lines) states. Recall that the observations from benchmark case B1 are at 1/4096 the resolution of the the simulation, since the velocity is observed at every eighth point in each spatial dimension and in time. The quality of the reconstruction is evident in the figure, with the predictions at the initial time capturing the large scales of the flow but appear to contain some small-scale, or high-wavenumber, deviations (panel $c$i). At the final time, in contrast, these small-scale deviations are mostly vanished and the reconstruction quality is, qualitatively, improved (panel $c$ii).

The convergence history of the cost function that generated the assimilated initial condition is shown in figure 4$a$. The normalization is performed using the cost function associated with advancing the initial guess which was obtained by spline-interpolation of the observations and their projection onto the divergence-free space. The monotonic decrease of the cost function is ensured by the accurate evaluation of its gradient using the discrete adjoint. After 100 L-BFGS iterations, the cost function is reduced to $2.8\%$ of its initial value.

A quantitative assessment of prediction accuracy starts with an evaluation of the root-mean-square errors between the predicted and true states,

where $\langle\bullet\rangle$ denotes averaging, and the subscript indicates the averaging domain, here over the volume $V$; the same convention is adopted throughout the work for errors $\mathcal{E}$ and correlation coefficients $\mathcal{C}$. The error is plotted in figure 4$b$ versus time. When errors and correlations coefficients between the predicted and true fields are reported, they are evaluated throughout the domain, and not only at observation locations. At the initial time, the spline-interpolation of the observations (red circle), projection onto the solenoidal space (cross) and the adjoint-variational state all have seemingly similar errors—the last being obtained after 100 L-BFGS iterations. A mild reduction in the initial errors is achieved by the projection, and a further modest improvement is achieved by the variational approach, but the initial errors remain approximately $4\%$ of the bulk velocity, which is of the same order of magnitude as the root-mean-square fluctuations. The important difference arises, however, when the initial conditions are advanced in time. When the interpolated observations are projected onto the divergence-free field and marched using the Navier-Stokes equations (dashed line), without any data assimilation, the errors amplify as expected due to the chaotic nature of the flow. At long times, after a transient divergence, non-linear effects become dominant and thus the estimation error saturates. Simply performing spatio-temporal interpolation of the observations would maintain lower errors, similar to the red circle, although would not be a solution to the governing equations. Now consider the errors when the initial condition is obtained from the adjoint-variational state estimation (solid line). The initial errors in the assimilated field decay with time, and the flow more closely tracks the trajectory of the true field in state space during the observation window (shaded region, $0\leq t\leq 4.5$). At $t=T$, the errors are an order of magnitude smaller than those from interpolating and advancing the initial condition or performing spatio-temporal interpolation of observations.

Beyond the observation window, the estimated state again diverges from the truth but remains more accurate than evolving the interpolated initial field during the interval $[T,3T]$. These results demonstrate the potential of the algorithm to provide a better prediction of the future state $t>T$, when observations are not available. If new data do become available at later times, the estimated state can be adopted as the initial guess, and the same variational procedure can be applied to drive the estimated state towards the true state, again. This idea is demonstrated in the figure: new observation data were provided in the interval $t\in[2T,3T]$, which is marked by the light green region. The adjoint variational algorithm was applied in that new interval, and the predicted state at $t=2T$ now yields a new trajectory that more closely follows the true flow (green line). While this point was noteworthy, the focus of the present effort is on characteristics of the state estimation in the first window $t\in[0,T]$, which are generally applicable to any subsequent interval of new observations.

The root-mean-square estimation error is evaluated in the horizontal plane, $\mathcal{E}_{xz}(q)$, are plotted as a function of wall-normal direction in figure 5. The errors in the interpolated initial guess (dashed lines) are proportional to the level of physical fluctuations in the velocity field and, as a result, the errors in the streamwise component are most dominant especially in the near-wall region where $u$-perturbations are most energetic. The estimated initial condition (thin black line) is slightly improved relative to the interpolated state. The key observation is, however, at $t=T$ where all three components of errors in the estimated state are an order of magnitude more accurate than advancing the interpolated state using the forward model. For the streamwise component, the error is lower than $0.5\%$ of the bulk velocity, or equivalently, $8\%$ of the peak value of root-mean-square streamwise fluctuation. Figure 5$b$ also shows the correlation coefficient,

At $t=T$, the estimated and the true fields are nearly perfectly correlated, which highlights the capacity of the assimilated field to reproduce the true trajectory of the flow in state space.

In order to explain the improvement in accuracy during the observation time-horizon, we evaluate the spectra of the estimation error,

where $\hat{\mathbf{u}}$ is the Fourier transform of $\mathbf{u}$ in the streamwise and spanwise directions. The spectra of the errors are reported in figure 6$a$ as a function of the horizontal $(k_{x},k_{z})$ wavenumbers; also shown in panel $b$ is the spectrum of the true velocity field. Errors appear largest in the low-wavenumber modes, but they should be viewed relative to the high energy content in these modes in figure 6$b$. Normalized by the spectral density in the true field, low wavenumbers are better reconstructed, since they are encoded in the sparse observations. A more important observation is the behaviour of the high-wavenumber components of the errors. The estimated initial condition (panel $a$i) has appreciable errors in those wavenumbers (higher than the interpolated state—not shown). However, since most of these modes decay with time (panel $a$ii), they have little impact on the estimation quality later within the assimilation time horizon. In addition, due to their rapid decay, these high-wavenumber initial errors do not appreciably affect the value of the cost function which is integrated over the entire observation window. As a result, they persist in the initial condition during the optimization procedure. An effective strategy to reduce these initial high-wavenumber errors is to incorporate time-dependent weights in the cost function, that amplify the contribution of early observations near $t=0$ (Wang et al., 2019a). One should caution, however, that not all the initial high-wavenumber errors are stable and decaying; Small components of that error are unstable and amplify at the Lyapunov rate, and although they are not perceptible within the present time horizon, they will dominate at longer times.

Consistent with the above integral and spectral measures of the estimation quality, the reconstructed velocity field at $t=T$ is almost identical to the true one (side views in figure 7). An accurate estimation of vortical structures is, however, more challenging since the computation involves velocity gradients. These structures are visualized using the $\lambda_{2}$ vortex identification criterion (Jeong & Hussain, 1995) and compared in figure 7 (gray isosurfaces). The prediction of vortical structures is very compelling, both in the near-wall and outer regions. In the former, the vortical tubes are attached to the wall and elongated in the streamwise direction; Farther from the wall, the lifted vortical tubes break down and generate small-scale structures that are mostly captured by the estimated state. In practice, reconstructing the vorticity field from coarse-resolution experimental data is a challenge (Suzuki, 2012). The satisfactory estimation quality in panel ($b$) demonstrates the potential of our algorithm to augment under-resolved turbulent data from experiments.

In the previous section, the observation data were free of any noise. In practice, however, experimental measurements invariably contain errors and, as such, may violate the governing equations, lead to statistical errors, and severely preclude accurate evaluation of derivatives especially in turbulent flows where strong velocity gradients are expected (see e.g. Bardet et al., 2010). For example, Abrahamson & Lonnes (1995) assessed the ability of conventional circulation and least-square methods to reproduce the vorticity field from randomly perturbed DNS velocity field. When $5\%$ noise was superposed onto the fully-resolved velocity data, the uncertainty in the computed vorticity field reached $40\%$, and most of the small-scale structures were absent in the reconstruction.

In order to evaluate the performance of adjoint-variational algorithm with noisy data, we contaminate the observed velocities by Gaussian random noise with standard deviation that is proportional to the local velocity component (5). The spatiotemporal resolution of the data and the estimation window $T^{+}=50$ remain the same as in the benchmark case. We consider three levels of noise: $\sigma=\{0,5,10\}\%$, which correspond to cases B1, N1 and N2 in table 2. The first guess of the initial condition is interpolation of of the noisy data at $t=0$.

After 100 L-BFGS iterations, the cost function is decreased to $\{2.8,27,45\}\%$ for the three levels of contamination. Since the cost function is defined as the difference between the estimated and contaminated data, and the latter do not satisfy the Navier-Stokes equations, the cost function cannot decrease to zero. The root-mean-square error of the estimated state relative to the true, uncontaminated flow field was evaluated and shows a similar decay from the initial to the final time as in the benchmark case without noise. Therefore, only the results at $t=T$ are examined here (figure 8). The estimation error increases with the noise level (from light gray to black lines), but remains within $2\%$ of the bulk velocity. Note that due to the mean flow, the observation noise in the streamwise direction can exceed $\sigma$ times the bulk velocity, which means that the estimation accuracy of the streamwise velocity actually exceeds the precision of observation data. Comparatively, the estimation errors of the wall-normal and spanwise velocity components are bounded by the observation noise, which is approximately $\sigma$ times the root-mean-square turbulence fluctuations. Overall, even with the highest noise level ($\sigma=10\%$), the correlation coefficient (12) between the estimated and true state is close to unity at all the $y$ locations, as shown by the red lines in figure 8$b$.

The reconstructed vortical structures are visualized in figure 9. Although some of the small-scale structures are not captured, most of the reconstructed wall-attached and detached vortical structures remain nearly identical to the true flow, and the estimation quality is almost independent of the noise level. A quantitative assessment of the quality of the vorticity field is provided in figure 10. When noisy observations are interpolated and vorticity is evaluated, the error (dashed lines) in the near-wall region reaches $10-40\%$ of the mean vorticity at the wall. The error of the vorticity field from adjoint-variational estimation (black solid lines) is within $4\%$ of the wall vorticity, and the estimation accuracy is robust against observation noise. These results demonstrate the superiority of adjoint-variational approach for evaluating velocity gradients and its robustness to observation noise.

Although the results thus far demonstrated the accuracy of the flow reconstructions, it is expected that the estimation quality depends on the spatiotemporal resolution of observations. And it is also of interest to query the lowest resolution requirement for accurate estimation. In homogeneous isotropic turbulence, it has been reported that reconstruction of the full field is successful only if the resolution of spatial observations satisfies $\Delta_{m}<15\eta$ (Li et al., 2020), where $\eta$ is the Kolmogorov lengthscale. An equivalent criterion has not, however, been proposed for anisotropic, wall-bounded turbulence where the effects of mean shear, advection and the no-slip boundary may alter the resolution requirements of observations. Hereafter, we revert to adopting noise-free data and focus on the influence of spatio-temporal resolution of observations on the accuracy of state estimation within the time horizon $T^{+}=50$.

We first consider the impact of spanwise spacing of observations that are fully resolved in the $x$-$y$ plane (cases RZ1$-$RZ4 in table 2). With the most coarse observations ($\Delta z_{m}^{+}=112$, case RZ4), the estimated state is visualized in figure 11 and compared to the true one. At the observation locations ($z=z_{m}$), velocity data are reproduced by the algorithm. At the midpoint between observation planes ($z=z_{m}+0.5\Delta z_{m}$), however, the estimation accuracy is notably compromised.

Since our interest is in the distribution of errors between observation planes, the estimation error is phase averaged in the span in addition to averaging in the streamwise direction, and is denoted $\mathcal{E}_{xz_{m}}(q)$. The results for cases RZ1$-$RZ4 are shown in figure 12. Only the error for $u$ component is plotted, and the results for $v,w$ components are similar. With the most poorly resolved data (panel $a$), the estimation error increases by an order of magnitude from the observation planes to the midpoint between them. The error in the near-wall region is of the same order of magnitude as the local turbulence fluctuations, which is significantly higher than the error at the channel center. As better-resolved observations are adopted (panels $b$-$d$), the inhomogeneity of the error distribution in the spanwise direction becomes weaker.

The effect of spanwise data resolution can be explained by examining the spanwise two-point correlations $\mathcal{R}_{u}(\Delta z)$ in figure 13($a$i). Since the dominant structures in the wall layer are streaks and streamwise vorticies that are narrow in the span, the two-point correlation decays faster at locations closer to the wall. As a result, the domain of dependence of the observation planes shrinks from channel center to the wall, which explains the high estimation error in the wall layer. Figure 13($a$ii) shows the reconstruction quality at $t=T$ in terms of the correlation coefficient between the true and estimated state from the least-resolved data ($\Delta z_{m}^{+}=112$). The profiles are qualitatively similar to the two-point correlations within $\Delta z^{+}<20$. However, while the turbulent structures decorrelate at larger distances, the accuracy of the reconstruction remains relatively higher, and $\mathcal{C}_{xz_{m}}(u^{\prime})$ returns to approximately unity as we approach the next observation plane at $z^{+}-z^{+}_{m}=112$.

Similar to the notion of the Taylor microscale,

we introduce a lengthscale for the domain of dependence of one observation plane by replacing $\mathcal{R}_{u}$ in (14) by the correlation coefficient $\mathcal{C}_{xz_{m}}(u^{\prime})$. The resulting lengthscale, which we denote $\Lambda_{Cz}$, is representative of the domain of accurate estimation. Both $\Lambda_{z,u}$ and $\Lambda_{Cz}$ are plotted in figure 13$b$, and have similar values across the height of the channel: The spanwise size of the domain of dependence of observations is similar in the Taylor microscale. Therefore, the criterion

must be satisfied to guarantee an accurate estimation of the local $u$ field. Similar criteria can be adopted for accurate estimation of $v$ and $w$ components using the Taylor microscales $\Lambda_{z,v}$ and $\Lambda_{z,w}$, respectively. Since those lengthscales are commensurate with $\Lambda_{z,u}$, the condition (15) suffices. Using (15), we can also interpret the estimation results with different spanwise data resolutions (figure 12). When $\Delta z^{+}_{m}=56$ (figure 12$b$), the criterion (15) is satisfied for the bulk of the channel ($2\Lambda^{+}_{z,u}\approx 80$, c.f. figure 13$b$) and starts to be violated for $y^{+}<70$. As such, while figure 12$b$ reports high prediction accuracy in the bulk ($\mathcal{C}_{xz_{m}}(u^{\prime})=0.97$), errors increase in that near-wall region ($\mathcal{C}_{xz_{m}}(u^{\prime})=0.83$) and become inhomogeneous in the span. When $\Delta z^{+}_{m}=28<\min_{y}2\Lambda_{z,u}^{+}$ (figure 12$c$), every point in the flow is covered by the domain of dependence of observations, so the estimation quality becomes more accurate and uniform at all the $y$ locations ($\mathcal{C}_{xz_{m}}(u^{\prime})=0.99$).

We compare the criterion (15) to the one for homogeneous isotropic turbulence (Yoshida et al., 2005; Lalescu et al., 2013; Li et al., 2020) that $\Delta_{m}\lesssim 15\eta$. The Kolmogorov scale in our case is $\eta\equiv(Re^{3}\mathcal{D})^{-1/4}$, where $\mathcal{D}=(2/Re)\langle s^{\prime}_{ij}s^{\prime}_{ij}\rangle$ is the viscous rate of dissipation and $s^{\prime}_{ij}=(\partial_{i}u_{j}^{\prime}+\partial_{j}u_{i}^{\prime})/2$ is the fluctuating strain rate tensor. The criterion $15\eta$ is plotted as red dashed line in figure 13$b$. Since the average dissipation is affected by the streamwise elongated structures in the channel, $\eta$ is larger than the spanwise size of the smallest eddies. Nonetheless, the trend is similar to the Taylor microscale condition provided above. Physically, both criteria demonstrate that the critical data resolution for accurate estimation of entire flow is within the transition zone between inertial and viscous dissipation ranges.

Next, we consider fully resolved observations in the cross-flow $y$-$z$ planes, and under-resolved in the streamwise direction (cases RX1$-$RX4 in table 2). Before examining the state estimation results, we report the Taylor microscales $\Lambda_{x,i}$ ($i=u,v,w$) in figure 13$c$, which is computed from the streamwise two-point correlations. Due to the near-wall elongated streaky structures, $\Lambda_{x,u}$ is largest among all three components, peaks near the wall and decays towards the channel center. We therefore expect that the estimation, particularly of the $u$ component, should be most accurate in the inner layer relative to the accuracy in the outer flow.

The root-mean squared estimation error is plotted in figure 14, where $\mathcal{E}_{x_{m}z}(u)$ is averaged in the span and phase-averaged in the streamwise direction, and also normalized by the r.m.s. fluctuations. Overall, the estimation error decreases as better-resolved data are included (panels $a$-$d$). For each observation resolution, the estimation quality deteriorates from the wall to channel center, as expected. Two notable differences from the effect of spanwise resolution are observed: (i) the estimation error is not symmetric with respect to the midpoint between observation planes, especially in panels ($a$-$b$); (ii) although the separation of observation planes in figure 14$a$ is larger than $2\Lambda^{+}_{x,u}(\approx 118)$ in the bulk (c.f. figure 13$c$), the estimation error does not increase appreciably between streamwise observation locations and remains within approximately $10\%$ of the local root-mean-squared fluctuations ($\mathcal{C}_{x_{m}z}(u^{\prime})\geq 0.95$ while $\mathcal{R}_{u}(\Delta x^{+}=96)\approx 0.6$). Both points are caused by the mean advection in the streamwise direction. Conceptually, every instant when data are recorded in the cross-flow $y$-$z$ plane corresponds to an accurately estimated “layer” in the spatiotemporal evolution of the flow, propagating by the advection velocity $U_{a}$. The thickness of such layers is approximately the Taylor microscale $\Lambda_{x,u}$ and the distance between two adjacent ones is approximately $U_{a}\Delta t_{m}$. Since the observation data are temporally well resolved, $U_{a}\Delta t_{m}\ll 2\Lambda_{x,u}$, the accurately estimated layers overlap with one another and lead to commensurate accuracy between observation locations. An example that further highlights this conceptual interpretation is considered next, where the temporal sampling rate is sufficiently low ($U_{a}\Delta t_{m}>2\Lambda_{x,u}$) in order to distinguish the accurately predicted layers associated with different ($x_{m},t_{m}$).

Consider the same coarse spatial resolution from case RX4 ($\Delta x^{+}_{m}=192$), and we now adopt a long time between observations, $\Delta t^{+}_{m}=7.4$. This case is denoted RT4, and a temporal evolution of the estimation error from $t=0$ to $t=\Delta t_{m}$ is shown in figure 15. Three layers of small errors are observed: ($A$) originates at the observation station at $t=0$ and advects downstream; ($B$) emerges between observation locations and times, such that the field accurately reproduces measurements at a downstream observation point at a subsequent measurement time $t=\Delta t_{m}$; ($B^{\prime}$) is similar to ($B$) but reaches the observation station at $t=2\Delta t_{m}$. Once ($A$) leaves the observation plane, the error increases, at least initially, due to the chaotic nature of turbulence and the absence of a nearby observation station to correct the field—behaviour later in time is discussed below. By comparison, the error in ($B$) decays with time as it approaches the observation location at time $t=\Delta t_{m}$.

A detailed space-time representation of the estimation error at the channel center is shown in figure 16. All the regions of small errors (blue) belong to one of the three classes described above, and can be associated with at least one observation station (black pluses). The pattern of the errors is repeated, anchored at $(\Delta x_{m},\Delta t_{m})$ and inclined according to the advection speed. Since $3U_{a}\Delta t_{m}\approx 2\Delta x_{m}$, regions ($A_{i}$) that originate from $(x,t)=(x_{m},0)$ will reach another observation station at $(x,t)=(x_{m}+2\Delta x_{m},3\Delta t_{m})$; the errors in that region thus initially increase moving away from $(x_{m},0)$ and reduce as that later observation point is approached. For a general spatiotemporal resolution of observations, regions ($A_{i}$) may not coincide with another observation position and time, and hence the errors would not undergo the later decay. From figure 16, it is evident that an accurate estimation of the entire flow state can be achieved by refining either the spatial or temporal data resolution such that the low-error regions overlap. This view is demonstrated in figure 17, where panels ($a$-$d$) correspond to increasing temporal resolution (cases RT4$-$1 in table 2). When the true state is observed sufficiently frequently, estimation error in the bulk region approaches homogeneity in the streamwise direction (panels $c$, $d$). Specifically, accurate reconstruction is achieved when the distance between the thin layers that sample observation points is within the domain of dependence of observations,

Recalling that $\Lambda_{x,v}\approx\Lambda_{x,w}<\Lambda_{x,u}$ (c.f. figure 13$c$), the equivalent criterion for accurate estimation of $v$ and $w$ is more restrictive. When the observations only satisfy the bound for $u$, the reconstruction is less accurate for the other two velocity components (e.g. $\mathcal{C}_{x_{m}z}(u^{\prime},v,w)=\{0.95,0.86,0.79\}$ when $\Delta x^{+}_{m}=192$ and $\Delta t^{+}_{m}=7.4$).

In summary, to guarantee an accurate estimation of local velocity component $i$, the streamwise and temporal data resolution must satisfy either one of two conditions:

The first condition is similar to the one for the spanwise data resolution, namely that the Taylor microscale must be resolved by observations. Note that it is supplemented by a statement regarding the observation time, which must be smaller than the Lyapunov timescale $\tau_{\sigma}$ because even a near perfect estimate of the flow state will diverge from the true solution over that period unless additional data are available for assimilation. The second condition is a reinterpretation that accounts for mean advection: Should the temporal resolution resolve the advected Taylor scale, the only requirement for spatial resolution becomes a condition based on the travel distance within the Lyapunov time. These two conditions are plotted in figure 18. We performed additional tests with $\Delta x_{m}^{+}=48$ and $\Delta t^{+}_{m}=\{0.45,0.90,1.80,3.70,7.40\}$. When the correlation coefficient between the true and estimated state was higher than 90%, the outcome was deemed successful; these cases are marked by circles. Inaccurate, or unsuccessful, reconstructions are marked by crosses. The outcomes of all the cases agree with the above two conditions.

In all the above cases, observations were distributed throughout the wall-normal extent of the channel, spanning both the inner and outer layers. In reality, however, experimental observations become progressively more difficult to obtain near the wall (Hutchins et al., 2005; Smits et al., 2011). And since most of the turbulence kinetic energy is produced near the wall (peak production at $y^{+}\approx 12$) and transported into the bulk region, lack of observation in this region may severely compromise the estimation of the full state. We therefore explore the accuracy of reconstructing the flow within the inner region when observations are only available in the outer layer and directly on the wall. Previous efforts for this configuration are limited and have generally focused on use of linear models (Baars et al., 2016; Illingworth et al., 2018). The present test case is the first attempt to perform the reconstruction based on the full nonlinear Navier-Stokes equations.

Specifically, we consider two types of observations: (i) sub-sampled velocity data available in the outer layer, $y^{+}>30$, with the same spatial and temporal resolution as the benchmark case; (ii) shear stresses $\boldsymbol{\tau}=(\tau_{xy},\tau_{zy})$ on both channel walls with the same resolution as velocity data. The observation time horizon is the same as the benchmark case, $T^{+}=50$. These observations are weighted in the cost function,

where $\|\cdot\|_{w}$ represent integration over top and bottom walls and the weight $\beta=1/5$. The choice of weighting parameter $\beta$ may be motivated by different objectives, such as balancing the two contributions in the cost function or their gradients or minimizing the condition number of the Hessian of $\mathcal{J}$. Here we adopt $\beta=1/5$ such that the stress term is commensurate with the velocity at the first fluid grid point that was observed in our benchmark case. The estimated state is compared with the benchmark case to highlight the impact of missing observations in the wall layer.

The predicted statistics without wall-layer observations, evaluated at $t=T$, are plotted in figure 19 and are compared with the true statistics. The shown quantities are horizontally averaged; the corresponding turbulence kinetic energy equation is,

The estimated profiles are accurate, even within the viscous sublayer and buffer layer which are void of observations and where the production and dissipation predominate. Although the interpolated initial condition also converges to the true statistics after a long time, our estimated state reaches the statistically stationary state faster, especially for dissipation. The errors in the estimated instantaneous velocity fields are shown in figure 20. Compared with the benchmark case (gray lines), the estimation quality in the outer layer is almost unaffected when the wall-layer data are removed. In the near-wall region, however, the estimation error reaches a maximum and exceeds twice the error in the benchmark case. The streamwise component is the most poorly reconstructed. Without data in the wall layer, the data assimilation algorithm starts with an interpolated velocity field as a first guess, which underestimates the mean flow near the wall; the final prediction of the algorithm retains a similar deficit in the mean, and hence the estimation error is largest in the streamwise component (black lines in figure 20). At $t=T$, the highest estimation error is approximately $4\%$ of the bulk velocity, which is smaller than the $16\%$ local root-mean square streamwise turbulence fluctuations. In addition, the correlation coefficient between the true and estimated state at $t=T$ is above 0.8 for all three components of the velocity field.

The wall layer is host to coherent structures, such as streaks and streamwise vorticies, that play an important role in the dynamics of the near-wall turbulence cycle (Hamilton et al., 1995) and in flow control. For example, the streak spacing must be resolved by a controller in order to effectively relaminarize turbulent channel flow (Sharma et al., 2011). Accurate prediction of these structures is therefore important. The error of our estimation at different lengthscales is reported in figure 21. The most accurately reconstructed structures are long in the streamwise direction ($\lambda_{x}\approx O(L_{x})$) and large waves in the span ($\lambda^{+}_{z}\in(10^{2},10^{3})$), which are typical features of the near-wall streaks and streamwise vorticies (Jiménez, 2018). The estimation error of these structures are mildly affected by removing the wall-layer data, due to their coherence across wall-normal locations, which facilitates the estimation even without near-wall observations. In contrast, reconstruction of small-scale near-wall structures is more sensitive to the lack of observations in that region, which is symptomatic of the weaker sensitivity of the outer flow to these structures.

Despite the overall increase of estimation errors when wall-layer data are not available, the instantaneous visualization of the coherent structures can still be compelling. Figure 22 shows a visualization of the predicted instantaneous streamwise fluctuation velocities for two predictions: without and with observations in the wall layer. The estimated streaky structures without wall-layer observations (panel $a$i) are irregularly spaced in the spanwise direction and meander downstream, as expected from a realistic channel flow. The zoomed in view (panel $a$ii) also includes the true state (line contours), which confirms that the state estimation is predictive: the true streaks are reproduced. The estimation accuracy is in fact comparable to the case when observations are available in the wall layer (panel $b$). The successful estimation of the missing wall layer is directly tied to the sensitivity of the observations in the outer flow and at the walls to the state in this region at the initial state. In other words, accurate reconstruction of the inner layer at the initial time is indispensable to match future observations far from the wall.

In practice, control of wall-bounded turbulence may rely on sensing and actuation at the wall, and effective strategies are often predicated on decoding the wall signature of outer turbulence structures. It has been demonstrated both mathematically (Lighthill, 1963; Constantin et al., 2011) and numerically (Eyink et al., 2020a, b) that all the interior vorticity is generated at the wall. The converse problem, specifically to what extent the initial flow state can be predicted from wall observations, has not been addressed as comprehensively. Previous efforts at $Re_{\tau}=100$, and using a variety of approaches, yielded an estimated state that is nearly uncorrelated with the true flow beyond the buffer layer (Bewley & Protas, 2004; Suzuki & Hasegawa, 2017; Liu & Hasegawa, 2020). The modest Reynolds number in those studies does not support the presence of large-scale outer structures that appreciably influence the wall stress fluctuations (Abe et al., 2004). Linear stochastic estimation (LSE) has been applied at higher $Re_{\tau}$ to reconstruct large-scale structures in the log layer from wall observations (Encinar & Jiménez, 2019). However, the predicted state from LSE does not satisfy the Navier Stokes equations and is not suitable for forecasting —a further discussion of LSE is provided below. Here we focus on adjoint-variational state estimation, specifically on the accuracy of the predicted state and the domain of dependence of wall observations. The dependence on Reynolds number will be examined for $Re_{\tau}=\{100,180,392,590\}$. The computational domain and grid resolution are summarized in table 3. Note that a smaller domain size is adopted for the $Re_{\tau}=\{392,590\}$ cases due to the limited computational resources.

The observations are fully-resolved shear stresses $\mathbf{\tau}=(\tau_{xy},\tau_{yz})$ and pressure $p$ at both walls, similar to those adopted by Bewley & Protas (2004) at the lower Reynolds number. The corresponding cost function is defined as,

For all the different $Re$ cases, the estimation window is the same in viscous units, $T^{+}=50$, which is long enough for perturbation in the bulk region to affect wall signals. The first guess of the initial condition is a linear stochastic estimation of the flow using observations at $t=0$. The stochastic estimator was constructed with a completely independent time series. All the estimated flows are obtained after one hundred L-BFGS iterations of the state-estimation 4DVar algorithm.

The correlation coefficients of the estimated and true fluctuating velocities at $t=T$ are shown in figure 23. The $Re_{\tau}=100$ results are comparable to those in figure 6$a$ from Bewley & Protas (2004) at the same Reynolds number. As $Re_{\tau}$ is increased, the correlation near the wall slightly deteriorates but remains sufficiently high to provide confidence in predictions. Precisely, for all Reynolds numbers, $C_{xz}>0.8$ when $y^{+}<15$. The correlations start to decay beyond the buffer layer, and those for $v$ and $w$ (panels $b$&$c$) do so nearly monotonically and with a similar slope for all $Re_{\tau}$. This trend highlight the challenge of interpreting turbulent flows from wall observations: the accurately predicted near-wall layer is a diminishingly smaller physical region as the Reynolds number is increased. An noteworthy exception is recorded in the correlation coefficients of the streamwise velocity (panel $a$) at $Re_{\tau}=\{392,590\}$: The initial decay outside the buffer layer is followed by a plateau within ($30<y^{+},y<0.3$) where the reconstruction remains marginally accurate because the large-scale structures in that region superimpose a footprint on wall signals (Abe et al., 2004; Mathis et al., 2009).

A comparison of the adjoint-variational and linear-stochastic estimated initial conditions is provided in figure 24. Both initial states are advanced using the Navier-Stokes equations, and their respective correlation coefficients with the true flow are reported as a function of time. At early times, both accurately reproduce the near-wall layer and outer large-scale $u^{\prime}$-structures, although the adjoint approach yields slightly better accuracy. At later times, the evolution of the LSE state diverges from the true flow, especially in the near-wall region. By contrast, the adjoint-variational estimation shadows the true trajectory in the near-wall layer and the large-scale motions in the outer layer.

Instantaneous views of the perturbation field for $Re_{\tau}=590$ are provided in figure 25, at $t=T$. In the outer layer, the estimated field matches the mean flow and large-scale motions, while the small-scale fluctuations in the bulk and vorticies that are detached from the wall are not captured. These results indicate that the wall signature at these Reynolds numbers is not sensitive to the wall-detached motions in the initial condition—a demonstration of the inherent difficulty of turbulence reconstruction from wall observations.

In light of the importance of the near-wall turbulence regeneration cycle, we turn to the detail of the true and estimated states in that region. The near-wall streaks (bottom $x-z$ planes in figure 25) are accurately reconstructed, and many of the vortical structures (gray isosurfaces in figure 26) are reproduced in the estimated field. Sheng et al. (2009) reported that violent near-wall ejections and sweeps contribute appreciably to the local turbulence energy. Examples of both events are shown in panels (i) and (ii), in figure 26: the true vortex lines above the stress local extrema are closely followed by the estimation, which demonstrates that the extreme events in the buffer layer are encoded in the wall signatures.

The accuracy of the estimated state can be viewed against the backdrop of the near-wall “autonomous cycle” (Jiménez & Pinelli, 1999; Jiménez & Moser, 2007). By artificially removing the outer flow in direct numerical simulations, Jiménez & Pinelli (1999) demonstrated that the near-wall dynamics, especially the regeneration cycle of streaks and vorticies, is self-sustained. And as observed in experiments, “ejection” and “sweep” events in that cycle have a wall signature in the form of minima and maxima of the wall shear stress (Sheng et al., 2009). Here we demonstrate that the wall stress is in fact encoded with the entire dynamics of the regeneration cycle, and our state estimation algorithm decodes these observations to discover the entire flow within the wall layer. The algorithm can also reconstruct the outer large-scale motions when their footprint is encoded in the wall measurements.