1. Introduction
Accurate predictions of turbulence must contend with its chaotic and multiscale nature: a slight deviation in the initial or boundary conditions exponentially amplifies and leads to an inaccurate prediction of all scales, and, hence, significant deviations of the trajectories in state space (Deissler Reference Deissler1986; Nikitin Reference Nikitin2018). In addition, it is difficult to precisely measure all the scales of turbulence, especially at high Reynolds numbers due to the larger separation of scales. Data assimilation methods aim to reconstruct the full state of a dynamical system from limited observations (Law, Stuart & Zygalakis Reference Law, Stuart and Zygalakis2015; Zaki & Wang Reference Zaki and Wang2021). These methods provide estimated trajectories that shadow the true state within an observation time horizon and that forecast the flow evolution beyond the observations, although, naturally, with progressively decreasing accuracy due to chaos. In homogeneous isotropic turbulence, recent studies have demonstrated that observing all the scales that are larger than approximately twenty Kolmogorov length scales is sufficient to recover all the missing smaller eddies, accurately – a phenomenon termed synchronization of chaos (Yoshida, Yamaguchi & Kaneda Reference Yoshida, Yamaguchi and Kaneda2005; Lalescu, Meneveau & Eyink Reference Lalescu, Meneveau and Eyink2013). The notion of synchronization is more restrictive than flow reconstruction: while the latter may refer to estimating a state that is close to the truth within a short time horizon, synchronization occurs only when the estimation error asymptotically approaches zero in time (Boccaletti et al. Reference Boccaletti, Kurths, Osipov, Valladares and Zhou2002). In the present work we investigate synchronization in turbulent channel flow, in physical rather than spectral space. Specifically, we report on the maximum size of a cloaked region that can be synchronized to the true flow trajectory using outer observations, and the required observations.
In order to incorporate observations into nonlinear dynamical systems, three types of data assimilation approaches have been developed: variational methods (Dimet & Talagrand Reference Dimet and Talagrand1986; Wang, Hasegawa & Zaki Reference Wang, Hasegawa and Zaki2019b; Li et al. Reference Li, Zhang, Dong and Abdullah2020), ensemble methods (Evensen Reference Evensen2009; Mons et al. Reference Mons, Chassaing, Gomez and Sagaut2016) and continuous data assimilation techniques (Charney, Halem & Jastrow Reference Charney, Halem and Jastrow1969; Yoshida et al. Reference Yoshida, Yamaguchi and Kaneda2005). Variational methods seek the optimal initial condition that minimizes a cost function, defined in terms of the distance between the estimated trajectory and available data. The minimization procedure requires the gradient of the cost function, which is efficiently computed using the adjoint equations. Ensemble methods do not involve an adjoint model. They consist of a prediction step where an ensemble of states are advanced using the forward equations, and an analysis step which assimilates observations to update the estimated flow state. The two classes can be combined, where an ensemble of forward simulations is adopted to construct a quadratic approximation of the cost function and to evaluate the gradient, which is generally referred to as the ensemble variational method (Liu, Xiao & Wang Reference Liu, Xiao and Wang2008; Mons, Wang & Zaki Reference Mons, Wang and Zaki2019). All of these techniques are widely adopted by the numerical weather prediction community and have been applied to estimating turbulent systems. However, due to the high-dimensional nature of turbulence, these methods are computationally expensive, require numerous, costly numerical simulations to obtain converged results (Jahanbakhshi & Zaki Reference Jahanbakhshi and Zaki2019; Wang, Wang & Zaki Reference Wang, Wang and Zaki2019a; Buchta & Zaki Reference Buchta and Zaki2021), which prohibits their application for synchronization of the estimated state with the true trajectory, especially at high Reynolds numbers.
Continuous data assimilation, which includes nudging techniques (Hoke & Anthes Reference Hoke and Anthes1976; Lakshmivarahan & Lewis Reference Lakshmivarahan and Lewis2013), augments the governing equations with a forcing term that drives the estimation towards available observations. Therefore, only one forward numerical simulation is required, which is the lowest cost compared with other approaches. Accurate predictions depend on the forcing scheme and, most importantly, the available observations. Construction of the forcing term is referred to as the ‘observer problem’ in control theory: given a dynamical system, the ‘observer’ refers to a forced system that assimilates limited data with the objective of synchronizing to the original system (Huijberts & Nijmeijer Reference Huijberts and Nijmeijer2001). Although different forcing strategies have been developed (Pogromsky & Nijmeijer Reference Pogromsky and Nijmeijer1998; Mohan, Liu & Vasudevan Reference Mohan, Liu and Vasudevan2017), they have a relatively minor influence on the success of synchronization in comparison to the amount of available observations of the true state. Therefore, we adopt the most straightforward and well-established approach: direct substitution where we replace part of the state vector by the corresponding observation data (Pecora & Carroll Reference Pecora and Carroll1990; Yoshida et al. Reference Yoshida, Yamaguchi and Kaneda2005). This approach can be viewed as a variant of the nudging technique in the limit of the forcing term having an infinitesimal relaxation time. The advantage of direct substitution against other nudging methods is the independence on any artificial parameters, and, thus, we can focus on the effect of observation on synchronization.
Recent efforts in homogeneous isotropic turbulence (Yoshida et al. Reference Yoshida, Yamaguchi and Kaneda2005; Lalescu et al. Reference Lalescu, Meneveau and Eyink2013; Vela-Martín Reference Vela-Martín2021) demonstrated that observations of the velocity field at all scales 
$k \eta < 0.2$, where 
$k$ is the wavenumber and 
$\eta$ is the Kolmogorov length scale, can guarantee synchronization of the smaller scales. This criterion is independent of the data assimilation techniques adopted, including direct substitution (Yoshida et al. Reference Yoshida, Yamaguchi and Kaneda2005), nudging (Clark Di Leoni, Mazzino & Biferale Reference Clark Di Leoni, Mazzino and Biferale2020) or adjoint-variational methods (Li et al. Reference Li, Zhang, Dong and Abdullah2020). Nikolaidis & Ioannou (Reference Nikolaidis and Ioannou2022) investigated synchronization in plane Couette flow by enforcing streamwise low wavenumber modes and obtained a similar criterion as in isotropic turbulence. No previous study has examined synchronization of wall turbulence in physical space or explored the anisotropy of the synchronization criterion along different flow directions. In channel flow the variation of the mean shear and streamwise advection all but eliminate the likelihood of a unified criterion across different locations and along different coordinates. We will therefore perform continuous data assimilation using direct substitution in physical space, and report on the maximum volume of turbulence that can be cloaked and yet synchronizes to the true trajectory by aid of observations of the rest of the system.
It is helpful to contrast synchronization in physical space to the previous studies that were performed in spectral space. There, the turbulence was prescribed at all wavenumbers above a cutoff value, and, hence, observations were available at any spatio-temporal location for all the energy-containing and inertial scales. Synchronization of a cloaked subvolume in physical space is markedly different because within the cloaked region, and for all times, none of the flow scales are known. Therefore, it is more difficult to anticipate whether the turbulence can be accurately generated from observations of the true state outside the subvolume. Consider, for example, removing a wall-attached horizontal layer in turbulent channel flow. This region includes the vorticity flux at the wall, which generates all the interior vorticity (Lighthill Reference Lighthill1963; Wu, Ma & Zhou Reference Wu, Ma and Zhou2007; Eyink, Gupta & Zaki Reference Eyink, Gupta and Zaki2020), and, therefore, synchronization using only outer observations may be difficult to achieve. If the layer includes the region of peak turbulence production, again whether it is possible to synchronize this ‘engine’ of turbulence is uncertain. Non-local interactions, for example, of outer large-scale motions with near-wall structures (Mathis, Hutchins & Marusic Reference Mathis, Hutchins and Marusic2009; Bernardini & Pirozzoli Reference Bernardini and Pirozzoli2011; Hwang et al. Reference Hwang, Lee, Sung and Zaki2016), may promote synchronization. Previous efforts have remarked on the role of these interactions in state estimate using linear (Baars, Hutchins & Marusic Reference Baars, Hutchins and Marusic2016; Illingworth, Monty & Marusic Reference Illingworth, Monty and Marusic2018) or nonlinear approaches (Sasaki et al. Reference Sasaki, Vinuesa, Cavalieri, Schlatter and Henningson2019; Wang & Zaki Reference Wang and Zaki2021), but none have examined synchronization. When the unknown layer is vertical and normal to the mean-flow direction, streamwise advection of upstream observations can aid synchronization. If the unknown vertical layer is, however, parallel to the flow, synchronization over the entire channel height is anticipated to be relatively more challenging. In addition, it is difficult to anticipate the amount of external observations that are required for synchronization to take place. All these considerations will herein be examined in the context of turbulent channel flow, for the first time.
In § 2 we provide details of the computational set-up, and introduce the direct substitution algorithm. We apply the method in § 3.1 to synchronize a horizontal wall-attached layer by observing the fully resolved outer flow, and we report the maximum layer thickness for successful synchronization and the influence on Reynolds number. Synchronization of wall-detached and of multiple layers are examined in §§ 3.2 and 3.4, respectively, and the equation for the evolution of synchronization error is derive in Appendix A. In § 3.5 we assess synchronization of a flow region using surface observations only, and discuss potential ways to improve the estimation accuracy of this approach. Additional tests for subdomain synchronization in Kolmogorov flow are presented in Appendix B. Concluding remarks are provided in § 4.
2. Computational approach
The reference flow configuration is a rectangular channel (see figure 1), which is periodic in the streamwise (
$x$) and spanwise (
$z$) directions, and bounded by two parallel no-slip surfaces in the vertical direction (
$y$). The channel half-height 
$h^*$ and bulk flow velocity 
$U_b^*$ are adopted as the reference length and velocity scales, where superscript 
$*$ denotes dimensional quantities. When quantities are scaled by inner variables, specifically the friction velocity 
$u_{\tau }^*$ and kinematic viscosity 
$\nu ^*$, they are marked by superscript 
$+$. The bulk and friction Reynolds numbers are 
$Re \equiv U_b^* h^*/\nu ^*$ and 
$Re_{\tau } \equiv u_{\tau }^*h^*/\nu ^*$, respectively.
Figure 1. Schematic of the reference and synchronization simulations of turbulent channel flow. A sample cloaked, or unobserved, horizontal region in the synchronization simulation is marked in green.
Our objective is to perfectly reconstruct the unknown turbulent field within a region 
$\varOmega _s$ of the channel, where we do not have any data, from observations of the outer flow in a region 
$\varOmega _f$. For example, in figure 1 the horizontal green layer 
$y \in [y_0,y_0+l_y]$ is the cloaked region 
$\varOmega _s$ and observations are available in the entire outer domain 
$\varOmega _f$. The full data assimilation domain is 
$\varOmega = \varOmega _s \cup \varOmega _f$. Perfect reconstruction of 
$\varOmega _s$ here means synchronization to the true flow that generated the outer observations, to within machine precision. We start by introducing the numerical schemes, computational set-up and synchronization algorithm.
The velocity field 
$\boldsymbol u$ and pressure 
$p$ satisfy the incompressible Navier–Stokes equations,
\begin{gather} \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{u} = 0, \end{gather}
\begin{gather}\frac{\partial \boldsymbol{u} }{\partial t} + \boldsymbol{\nabla} \boldsymbol{\cdot} (\boldsymbol{uu}) ={-} \boldsymbol{\nabla} p + \frac{1}{Re} \nabla^2 \boldsymbol{u}, \end{gather}
where 
$t$ represents time. These equations are spatially discretized on a staggered grid with a local volume-flux formulation (Rosenfeld, Kwak & Vinokur Reference Rosenfeld, Kwak and Vinokur1991), and advanced using a second-order fractional-step method (Kim & Moin Reference Kim and Moin1985; Bell, Colella & Glaz Reference Bell, Colella and Glaz1989) that adopts Adams–Bashforth scheme for the advection terms and Crank–Nicolson for diffusion. Due to periodicity in 
$x$ and 
$z$, the pressure Poisson equation is solved using Fourier transform in these directions and tri-diagonal inversion in the 
$y$ direction. The algorithm has been validated and applied in a number of direct numerical simulations of transitional and turbulent flows (Zaki Reference Zaki2013; Marxen & Zaki Reference Marxen and Zaki2019).
The true, or reference, system 
$\boldsymbol u_r(t)$ is an equilibrium turbulent flow that is sustained by a constant pressure gradient in the streamwise direction. The domain sizes and grid resolutions for all the considered Reynolds numbers are summarized in table 1. The two highlighted configurations, at 
$Re_{\tau }=\{590, 1000\}$, will be the focus of the majority of the discussion. The influence of domain size is examined at 
$Re_{\tau }=590$, where a larger domain 
$(L_x,L_z) = (8{\rm \pi}, 3{\rm \pi} )$ is also considered in order to accommodate very-large-scale structures (Abe, Kawamura & Choi Reference Abe, Kawamura and Choi2004; Del Álamo et al. Reference Del Álamo, Jiménez, Zandonade and D Moser2004). The time step size 
$\Delta t$ in each case is chosen to guarantee that the Courant–Friedrichs–Lewy number is less than one half.
Table 1. Computational domains and grid sizes for simulations at different Reynolds numbers.
The observer system, or synchronization simulation, is denoted 
$\boldsymbol u_s(t)$ and is performed in the domain 
$\varOmega = \varOmega _s \cup \varOmega _f$ which, unless otherwise stated, is the same as the reference simulation. Given the fully resolved true velocity field 
$\boldsymbol u_r(t)$ in subvolume 
$\varOmega _f$, these data are directly enforced onto the synchronization simulation, 
$\boldsymbol {u}_s=\boldsymbol {u}_r\ \forall \boldsymbol {x}\in \varOmega _f$ and 
$t$. The reference pressure field 
$p_r(t)$ is not observed, and, hence, the pressure in the synchronization simulation 
$p_s(t)$ differs from the true state throughout the entire domain 
$\varOmega$. Infusing the velocity observations into the synchronization simulation is performed using the direct substitution algorithm 1, which is also illustrated in figure 1. Due to the non-local nature of the Navier–Stokes solution, errors within 
$\varOmega _s$ instantaneously contaminate the solution within 
$\varOmega _f$, and, therefore, the direct substitution procedure is enforced every time step. Once applied, 
$\boldsymbol u_s$ differs from 
$\boldsymbol u_r$ only within the unobserved region 
$\varOmega _s$, where the synchronization error is defined,
\begin{equation} \mathcal{E}_{s} \equiv \langle \|\boldsymbol u_s - \boldsymbol u_r \|^2 \rangle_{\varOmega_s}^{1/2}, \end{equation}
where 
$\| \boldsymbol a\|$ is the 2-norm of vector 
$\boldsymbol a$, and 
$\langle \bullet \rangle _{\varOmega _s}$ denotes averaging over 
$\varOmega _s$. Theoretically, synchronization refers to
\begin{equation} \lim_{t \to \infty} \langle \|\boldsymbol u_s - \boldsymbol u_r \|^2 \rangle_{\varOmega_s}^{1/2} = 0 \end{equation}
for any arbitrary choice of initial condition 
$\boldsymbol u_s(0)$ within 
$\varOmega _s$. The condition (2.4) also guarantees that the pressure field 
$p_s$, which is not observed, converges to the reference state 
$p_r$ throughout the entire domain. Due to finite numerical precision and integration time, the above condition will be approximated by
\begin{equation} \mathcal{E}_s(t) < \beta, \quad \text{when} \ t < T. \end{equation}
The threshold 
$\beta$ is set to 
$10^{-15}$ due to our double-precision arithmetic, and 
$T=80$ is sufficiently long to achieve a conclusive trend.
The initial condition 
$\boldsymbol u_s(t=0)$ for 
$\boldsymbol {x}\in \varOmega _s$ is estimated using one of three approaches: (i) a trivial initial condition 
$\boldsymbol {u}_s(t=0) = \boldsymbol {0}$; (ii) the local mean velocity superposed with white noise proportional to the local root-mean-square (r.m.s.) fluctuations, 
$\boldsymbol u_s = \langle \boldsymbol {u}_r\rangle + \boldsymbol {u}_{\epsilon }$, where 
$u_{i,\epsilon } \sim N(0,u_{i,rms})$ is a Gaussian random field; (iii) a slight perturbation to the true state, 
$\boldsymbol u_s = \boldsymbol u_r + \epsilon \boldsymbol {u}_{\epsilon }$, where 
$\epsilon = 10^{-4}$. When synchronization occurs, all these initial conditions lead to an exponentially decreasing estimation error with the same rate, which is equal to the leading Lyapunov exponent of the sub-dynamics within 
$\varOmega _s$. If 
$\varOmega _s$ exceeds a certain size or if the available observations within 
$\varOmega _f$ are insufficient, synchronization is not possible. In such cases, even simulations starting from the initial estimate (iii), which is infinitesimally close to the true state, generate exponentially amplifying errors relative to the reference trajectory. Simulations of type (iii) are initially infinitesimally close to the true state, and, hence, are designed to probe the linear exponential stability of the subsystem; these computations will be termed Lyapunov experiments.
Algorithm 1: Direct substitution algorithm.
3. Results
3.1. Synchronization of a horizontal wall-attached layer
We start by attempting to synchronize a wall layer, 
$\varOmega _s = \{ \boldsymbol x \in [0,L_x] \times [0,l_y] \times [0,L_z]\}$, similar to the set-up in figure 1 with the lower boundary on the wall. In the limit 
$l_y\rightarrow 0$, synchronization is essentially trivial. As 
$l_y$ increases, successful synchronization becomes less inevitable: the absence of the true vorticity source at the wall from the observer simulation may impede synchronization, while the coupling of the outer observations within 
$\varOmega _f$ to the near-wall region may aid in the convergence of the simulation to the reference trajectory. Ultimately the balance of restoring and destabilizing effects is anticipated to lead to divergence of the sub-dynamics from the reference simulation, which in the limit of large 
$l_y$ is again not surprising since the entire streamwise and spanwise dimensions are eliminated, and all scales/wavenumbers in those dimensions at a given height are unknown. The interesting aspect of this problem is therefore not the limiting behaviours, but whether there exists a specific, or critical, value of 
$l_y$ across which the behaviour changes and how it depends on the Reynolds number.
A qualitative account of synchronization at 
$Re_{\tau }=1000$ is provided in figure 2, when 
$l_y^+ = 28$. The figure shows a side view of the channel, with the top panels displaying the reference simulation and the bottom panels focused only on the cloaked region 
$y \in \varOmega _s$ of the synchronization simulation. Panel (aii) is the initial estimate of 
$\boldsymbol {u}_s(t=0)$ within 
$\varOmega _s$, which is the mean flow plus white noise. Synchronization of the wall layer proceeds in two stages: during an initial transient (panel bii), the white noise decays, and the velocity field in the vicinity of the outer observations becomes qualitatively more accurate. However, quantitative accuracy is not yet achieved especially in the near-wall region. In the limit of long time (panel cii), the entire wall layer is accurately reconstructed. The synchronization process is examined in Fourier space in figure 3, where the pre-multiplied spectra of the streamwise velocity is reported. The colour and line contours correspond to 
$\hat {u}_s$ and 
$\hat {u}_r$, and 
$\hat {\bullet }$ denotes Fourier transform in the reported direction, e.g. streamwise direction in figure 3. Here too the two stages of synchronization are evident: immediately after the initial time (panel a), the small-scale fluctuations are overestimated due to the memory of white noise in the initial condition. As the outer observations are imposed at all scales, the estimated spectra converge to the truth at the same rate across all wavelengths (panels b,c).
Figure 2. Synchronization of a horizontal, wall-attached layer with 
$l_y^+ = 28$, at 
$Re_{\tau }=1000$. Contours are the instantaneous streamwise velocity fluctuations, calculated by subtracting the true mean, at 
$z = L_z/2$. Results are shown for (a–c) 
$t^+=\{ 0, 40, 160\}$; (i) true state; (ii) synchronization simulation. Dashed line: 
$y^+ = 28$.
Figure 3. Streamwise pre-multiplied spectra of the streamwise velocity, averaged in the spanwise direction, 
$\log _{10} \langle k_x^+ |\hat u^+|^2 \rangle _z$. Colours, synchronization simulation of a wall-attached horizontal layer with 
$l_y^+ = 28$; lines, reference simulation. (a–c) 
$t^+ = \{13,40,160\}$.
To quantify the accuracy of the synchronization simulations, we introduce the horizontally averaged error
\begin{equation} \mathcal{E}_{xz}(q) = \langle \left( q_s - q_{r} \right)^2 \rangle_{xz}^{1/2}, \quad q = u, v, w, p. \end{equation}
At 
$t=0$, the error 
$\mathcal {E}_{xz}(q)$ is approximately 
$\sqrt {2}$ of the r.m.s. fluctuations for all three components, because the superposed white noise in 
$\boldsymbol u_s(0)$ is uncorrelated with the true fluctuations. The temporal evolution of 
$\mathcal {E}_{xz}$ is plotted in figure 4, normalized by the local r.m.s. fluctuations. Unlike the velocity (panels a–c) where the errors vanish in 
$\varOmega _f$, the pressure is not observed and hence has finite error for 
$y^+ > 28$ (panel d). The two stages of error decay are evidenced in figure 4. At the initial stage 
$t^+ \sim O(10)$, the error remains highest near 
$y^+ = 0$ due to the delayed impact by observations. Beyond 
$t^+ = 100$, the synchronization errors within the entire removed layer are diminishing for both velocity and pressure, and uniformly tend to zero. As such, it is reasonable to focus on the volume-averaged error 
$\mathcal {E}_s$ (solid line in figure 5a), which decays until it reaches 
$10^{-15}$ dictated by our double-precision implementation.
Figure 4. Time and wall-normal dependence of the synchronization error. The error is averaged in the horizontal plane and normalized by the local r.m.s. fluctuation, 
$\log _{10} (\mathcal {E}_{xz}(q)/q^{\prime }_{rms} )$. Results are shown for (a–d) 
$q=\{u, v, w, p\}$; (a) 
$\log _{10} (\mathcal {E}_{xz}(u)/u^{\prime }_{rms} )$, (b) 
$\log _{10} (\mathcal {E}_{xz}(v)/v^{\prime }_{rms} )$, (c) 
$\log _{10} (\mathcal {E}_{xz}(w)/w^{\prime }_{rms} )$, (d) 
$\log _{10} (\mathcal {E}_{xz}(p)/p^{\prime }_{rms} )$.
Figure 5. (a) Temporal dependence of the volume-averaged synchronization error 
$\mathcal {E}_s$ normalized by the initial value 
$\mathcal {E}_{s,0}$. Results are for 
$Re_{\tau }=1000$: (dashed line) 
$l_y^+ = 8$; (dashed-dotted line) 
$l_y^+=18$; (solid line) 
$l_y^+=28$. (b) Effect of observation noise level on synchronization error averaged over 
$\varOmega _s$, when 
$l_y^+ = 28$. Black, blue, green: 
$\epsilon = \{0, 0.1, 0.5\}\%$.
We have also examined the effect of observation noise on synchronization. Specifically, we have contaminated the observed velocities which are imposed within 
$\varOmega _f$ by Gaussian random noise, proportional to the local r.m.s. fluctuations, 
$\boldsymbol {u}_s = \boldsymbol {u}_r + \epsilon \boldsymbol {u}_{\epsilon }\ \forall \boldsymbol {x}\in \varOmega _f$ and 
$t$, and where 
$u_{i,\epsilon } \sim N(0,u_{i,rms})$. For two levels of observation noise, 
$\epsilon =\{0.1, 0.5\}\%$, the synchronization errors within 
$\varOmega _s$ were evaluated and are reported in figure 5(b). The errors do not reduce to machine precision (blue and green curves), but rather saturate at a similar order of magnitude as the observation noise. It is most important to note that the synchronization rate remains unaffected by the noise, which demonstrates the robustness of the phenomenon.
As the thickness of the unknown layer 
$l_y$ is reduced, synchronization is precipitated at a faster rate, as shown in figure 5(a). The rate is evaluated from 
$\mathcal {E}_s(t)$ using an exponential regression
\begin{equation} \mathcal{E}_s = A \exp(\alpha t) = A \exp(\alpha^+ t^+), \end{equation}
where we have introduced the synchronization exponent 
$\alpha$. Only data within the range 
$10^{-1} \le \mathcal {E}_{s}/ \mathcal {E}_{s,0} \le 10^{-6}$ are used for the regression (3.2) in order to eliminate the transient effects at the beginning of simulations. The dependence of 
$\alpha$ on 
$l_y$ and on the Reynolds number is summarized in figure 6. The results demonstrate that indeed as 
$l_y$ increases, the synchronization rate reduces and ultimately becomes positive beyond a critical value. Note that, for the diverging 
$\alpha > 0$ cases, the exponent was determined using the Lyapunov-type experiments where the initial estimate 
$\boldsymbol {u}_s(t=0)$ was infinitesimally close to the true solution within the cloaked region 
$\varOmega _s$ (case iii, as explained in § 2).
Figure 6. Dependence of the synchronization exponents 
$\alpha ^+$ on the thickness 
$l_y^+$ of the cloaked wall-attached layer, at different Reynolds numbers. Green, 
$Re_{\tau }=180$; blue, 
$Re_{\tau }=392$; red, 
$Re_{\tau }=590$; black circles, 
$Re_{\tau } = 1000$. The dotted line marks 
$\alpha = 0$.
When scaled with viscous units, the profile of 
$\alpha ^+(l_y^+)$ is essentially independent of the Reynolds number up to 
$Re_{\tau } = 1000$. The critical value for 
$l_y$ that corresponds to zero growth rate,
\begin{equation} \alpha(l_{y,c}) = 0, \end{equation}
is 
$l_{y,c}^+ \approx 32$, which is the maximum thickness of the wall layer that can be removed without disrupting successful synchronization.
The quantitative relation in figure 6 is robust against the initial condition of the reference simulation 
$\boldsymbol u_r(0)$, the initial estimate of 
$\boldsymbol u_s(0)$ within the cloaked region 
$\varOmega _s$, and the global domain size. In particular, when the global domain size is increased from 
$(L_x,L_z)=(2{\rm \pi},{\rm \pi} )$ to 
$(L_x,L_z)=(8{\rm \pi},3{\rm \pi} )$ at 
$Re_{\tau }=590$, the value of the synchronization exponent for 
$l_y^+ = 18$ remains unchanged to within less than 
$1\,\%$. Therefore, synchronization of the wall layer is unaffected by the associated change in the very-large-scale structures in the outer flow (Marusic & Perry Reference Marusic and Perry1995; Marusic & Monty Reference Marusic and Monty2019).
The critical thickness 
$l_{y,c}^+ \approx 32$ demonstrates that the dynamics within the viscous sublayer and buffer layer are interpretable from outer observations. Even the instantaneous flux of vorticity at the wall is accurately reconstructed by enforcing the fully resolved outer flow field. Our results also provide a new perspective on the influence of outer large-scale structures on near-wall turbulence (Del Álamo et al. Reference Del Álamo, Jiménez, Zandonade and D Moser2004; Mathis et al. Reference Mathis, Hutchins and Marusic2009). Previous studies by Jiménez & Pinelli (Reference Jiménez and Pinelli1999) indicate that the outer velocity field near 
$y^+ = 60$ can sustain the near-wall cycle of streaks and streamwise vortical structures, but there was no guarantee of synchronization since they did not impose the velocity from a reference simulation. At that height, the predicted near-wall cycle adopts a different trajectory than when the entire channel is simulated. Our results further demonstrate that providing reference data and attempting to synchronize the near-wall layer is only possible for a thinner region, 
$l^+_{y,c} \approx 32$. Finally, we remark that 
$y^+ \leq 32$ is a diminishing small physical region as 
$Re$ is increased, which is indicative of an increasing difficulty of synchronization. The scales of turbulence that are commensurate with the critical thickness are discussed in the next section, where we examine the impact of placing the cloaked horizontal layer at different wall-normal heights.
3.2. Synchronization of a wall-detached layer
When the cloaked layer 
$\varOmega _s = [0,L_x]\times [y_0,y_0+l_y]\times [0,L_z]$ is detached from the wall, 
$y_0 > 0$, we can anticipate that the critical thickness for synchronization will change just as the scales of turbulence do, in particular the vertical size of the structures since 
$\varOmega _s$ will continue to exclude all the horizontal scales. In these tests, while the true vorticity flux at the wall becomes part of the observed region 
$\varOmega _f$, every new value of 
$y_0$ is associated with a shift in the dominant balance of turbulence dissipation, production and transport. In addition, the scaling of turbulence shifts from viscous units for 
$y^+ < 100$ (Lee & Moser Reference Lee and Moser2015) to inertia dominated for the large scales in the outer part (Adrian Reference Adrian2007; Smits, McKeon & Marusic Reference Smits, McKeon and Marusic2011). In order to determine the impact of 
$y_0$ on synchronization, we focus on two Reynolds numbers that provide sufficiently extended log regions, 
$Re_{\tau }=\{590, 1000\}$. The lower boundary of 
$\varOmega _s$ is gradually increased from 
$y_0^+ = 0$ up to 
$y_0^+ = 100$, which corresponds to removing part of the outer layer, 
$y_0=\{0.17, 0.10\}$ for 
$Re_{\tau }=\{590, 1000\}$, respectively. Synchronization of a horizontal layer located at the channel centre will also be briefly examined.
The profile of the synchronization exponent 
$\alpha ^+(l_y^+)$ at 
$y_0^+=100$ is shown in figure 7(a) (dashed line and crosses), and is compared with the exponent when 
$y_0^+=0$ (solid line and circles). When the cloaked layer is thin, 
$l_y^+ \lesssim 10$, the synchronization process in the log layer is slower than within the wall layer. As the removed region in the log layer is expanded, 
$\alpha ^+$ monotonically increases but at a relatively weak slope. As a result, a much thicker layer 
$l_{y,c}^+\approx 64$ is guaranteed to synchronize when 
$y_0^+ = 100$. The synchronization exponent 
$\alpha ^+(l_y^+)$ and the critical thickness 
$l_{y,c}^+$ in viscous units are unaffected by the Reynolds number (figure 7a), which is consistent with the inner scaling of the near-wall kinetic energy budget. This trend suggests that synchronization of a wall-detached layer up to 
$y_0^+ = 100$ is still governed by similar dynamical considerations.
Figure 7. (a) Synchronization exponent with (solid line, circle) 
$y_0^+ = 0$ and (dashed line, cross) 
$y_0^+=100$. Red lines, 
$Re_{\tau }=590$; black symbols, 
$Re_{\tau } = 1000$; horizontal dotted line, 
$\alpha ^+ = 0$. (b) Production (blue) and dissipation (green) of synchronization errors when (solid, dashed) 
$y_0^+ = \{0,100\}$. Only the results at 
$y_0^+ = 100$ are shown in the inset. (c) Symbols are the critical thicknesses 
$l_{y,c}^+$ as a function of the wall-normal height of the cloaked layer. Grey line: critical length scale in isotropic turbulence, 
$16 \eta ^+$. Blue, green, black: averaged Taylor microscales based on 
$\{u, v, w\}$.
Starting from the equations for 
$\boldsymbol u_r$ and 
$\boldsymbol u_s$, we derive the governing equation for the squared, volume-averaged synchronization error 
$\mathcal {E}^2_{s} = \langle \| \boldsymbol u_s - \boldsymbol u_r \|^2 \rangle _{\varOmega _s} : = \langle \| \boldsymbol e \|^2 \rangle _{\varOmega _s}$:
\begin{equation} \frac{1}{2} \frac{{\rm d}}{{\rm d} t}\langle \| \boldsymbol e\|^2\rangle_{\varOmega_s} ={-}\langle\boldsymbol{e}, \boldsymbol{e} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{r}\rangle_{\varOmega_s} -\nu \langle \| \boldsymbol{\nabla} \boldsymbol e\|^2\rangle_{\varOmega_s}+B, \end{equation}
where the notation 
$\langle \boldsymbol a, \boldsymbol b \rangle _{\varOmega _s} = \langle a_i b_i\rangle _{\varOmega _s}$ is the inner product between the two vectors, volume averaged over 
$\varOmega _s$. The first two terms on the right-hand side of (3.4) quantify the production and viscous dissipation of synchronization errors. The last term 
$B$ is due to the non-zero divergence of 
$\boldsymbol e$ on 
$\partial \varOmega _s$, i.e. on the boundaries of the unobserved region 
$\varOmega _s$. Since 
$B$ only contains surface integration (see Appendix A for details), it is generally much smaller than the production and dissipation terms. Similar equations have been derived by Henshaw, Kreiss & Yström (Reference Henshaw, Kreiss and Yström2003) for isotropic turbulence and adopted to estimate the critical cutoff wavenumber. Normalizing (3.4) by 
$\mathcal {E}_s^2$ and assuming 
$\mathcal {E}_s = A \exp {(\alpha t)}$, the equation for exponent 
$\alpha$ is
\begin{equation} \alpha = \frac{1}{\langle \| \boldsymbol e\|^2\rangle_{\varOmega_s}} (-\langle\boldsymbol{e}, \boldsymbol{e} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{r}\rangle_{\varOmega_s} - \nu \langle \| \boldsymbol{\nabla} \boldsymbol e\|^2\rangle_{\varOmega_s} + B ) \equiv \mathcal{P}_s - \mathcal{D}_s + \mathcal{B}_s, \end{equation}
where 
$\mathcal {P}_s$, 
$\mathcal {D}_s$, 
$\mathcal {B}_s$ are the normalized production, dissipation and boundary terms, respectively.
The qualitative behaviour of the synchronization exponent and critical thickness in figure 7(a,b) can be understood with reference to the production and dissipation of errors in (3.5). These terms are averaged in time after the initial transient, 
$\mathcal {P}_{st} = \langle \mathcal {P}_s \rangle _t$ and 
$\mathcal {D}_{st} = \langle \mathcal {D}_s \rangle _t$, and the results are reported in figure 7(b). For a thin layer (
$l_y \rightarrow 0$), viscous dissipation of the errors (green lines) dominates the balance. As a thicker layer of the flow is cloaked, the normalized dissipation significantly reduces and the normalized production rate of errors (blue) increases and exceeds dissipation beyond the critical thickness 
$l_{y,c}^+$. Despite a more pronounced dissipation of errors when the cloaked layer is wall attached (
$y_0^+ = 0$, solid green) compared with wall detached (
$y_0^+ = 100$, dashed green), the critical thickness is smaller due to the faster increase in production near the wall (solid vs dashed blue curves). The later effect is primarily due to the stronger mean shear, where 
$\partial \langle u_r\rangle / \partial y$ at 
$y_0^+ = 0$ is approximately forty times larger than the value at 
$y_0^+ = 100$.
The dependence of the critical thickness 
$l_{y,c}$ on distance from the wall is reported in figure 7(c). At each 
$y_0$, a bisection approach is adopted to identify the interval 
$[l_{y,n},l_{y,n} + \Delta y]$ (red bars in figure 7c) where the exponent 
$\alpha$ changes from negative to positive. The critical thickness (red dot) thus lies within this interval, and the value is estimated by linear interpolation using the exponents at 
$l_{y,n}$ and 
$l_{y,n} + \Delta y$. When the lower boundary of the cloaked layer is within the range 
$0 \le y_0^+ \le 12$, i.e. up to the height of peak kinetic energy production, the critical thickness for synchronization is practically unchanged. In the context of (3.4), this behaviour is due to a comparable balance between dissipation and production of errors when the layer starts within 
$y_0^+ \le 12$. However, as 
$y_0$ is more distant from the wall, a thicker layer of the flow can be synchronized to the reference state, or in other words the dynamics within 
$\varOmega _s$ can be discovered from outer observations.
The critical thickness 
$l_{y,c}$ can be compared with the synchronization criterion from isotropic turbulence, which was originally expressed in terms of a critical wavenumber 
$k_c = 0.2 / \eta$, where 
$\eta$ is the Kolmogorov scale. The corresponding resolution of observations in physical space is 
$\Delta _c \equiv {\rm \pi}/ k_c \approx 16\eta$. For channel flow, we evaluate 
$\eta = (Re^3 \mathcal {D}_r)^{-1/4}$, where 
$\mathcal {D}_r$ is the dissipation rate of turbulent kinetic energy, and the curve 
$16\eta$ is compared with 
$l_{y,c}$ in figure 7(c). Despite the qualitative similarity, the Kolmogorov length scale is not ideal for this comparison, and cannot account for the anisotropy of the flow when we consider cloaked layers with different orientations. A more suitable choice is the Taylor microscale which quantifies the size of flow structures along different directions. For synchronization in a horizontal layer, the wall-normal Taylor microscale 
$\varLambda _{y,i}$ is the most relevant:
\begin{equation} \varLambda_{y,i} = \left(\left.-\frac 12 \frac{{\rm d}^2 \mathcal{R}_i}{{\rm d} (\Delta y)^2} \right|_{\Delta y=0} \right)^{{-}1/2} \quad \text{for}\ i=u, v, w, \end{equation}
where 
$\mathcal {R}_i(\Delta y) = \langle u_i(y)u_i(y+\Delta y) \rangle _{xzt}$ is the wall-normal two-point correlation, averaged over the homogeneous horizontal directions and time. The Taylor microscale 
$\varLambda _{y,i}$ quantifies the wall-normal size of flow structures, and has recently been related to the domain of sensitivity of velocity observations (Wang & Zaki Reference Wang and Zaki2021). Therefore, we can use the Taylor microscales at 
$y_0$ and 
$y_0 + l_y$,
\begin{equation} 2 \overline{\varLambda_{y,i}} = \varLambda_{y,i}(y_0) + \varLambda_{y,i}(y_0+l_y), \end{equation}
as an estimate of the thickness of the cloaked layer that is entirely within the domain of sensitivity of the outer observations. The agreement between 
$2\overline {\varLambda _{y,w}^+}$ and 
$l_{y,c}^+$ in figure 7(c) suggests that the critical thickness is indeed proportional to the domain of sensitivity of the available turbulence data. A physical interpretation of the dependence on the Taylor microscale is provided, with reference to its definition in isotropic turbulence,
\begin{equation} \varLambda = \sqrt{15 \frac{\nu}{\mathcal{ D}_r}}\,u^{\prime}_{rms}. \end{equation}
The square-root term in (3.8) is proportional to the Kolmogorov time scale 
$\tau _{\eta } = \sqrt {\nu /\mathcal {D}_r}$. Therefore, physically, the Taylor microscale is a measure of the distance swept by Kolmogorov eddies during their lifetime, while advected by the root-mean-squared velocity fluctuations. As long as the flow data down to the Kolmogorov eddies can be swept from 
$\partial \varOmega _s$ throughout the cloaked region, prior to their dissipation, synchronization is guaranteed. In turbulent channel flow the lifetime of Kolmogorov eddies and the swept height increase with distance from the wall and, as a result, a thicker wall-normal layer can be decoded from outer observations.
The previous arguments remain applicable in the bulk region, which is evidenced by additional synchronization tests for a cloaked horizontal layer located at the channel centre, 
$\varOmega _s = \{\boldsymbol x \in [0,L_x] \times [1 - l_y/2,1+l_y/2] \times [0,L_z]\}$. The identified critical thickness is 
$l_{y,c}^+ \approx 140$, for 
$Re_{\tau }=590$, which is comparable to twice the Taylor microscale 
$2(\varLambda _{y,u},\varLambda _{y,v},\varLambda _{y,w}) \approx (139,163,110)$ at the boundary 
$y = 1 - l_{y,c}/2$. It is noteworthy that 
$l_{y,c}^+ \approx 140$ is of the same order of magnitude as the previously reported criterion for synchronization in homogeneous isotropic turbulence, 
$\Delta _c^+ \approx 16 \eta ^+ \approx 78$ at 
$y = 1 - l_{y,c}/2$ (Yoshida et al. Reference Yoshida, Yamaguchi and Kaneda2005). Despite this similarity, it is important to recall that the interpretation of the two criteria are fundamentally different. The results for isotropic turbulence are in terms of a critical wavelength in Fourier space. Specifically, the small-scale turbulence below 
$\Delta _c^+$ can be accurately reconstructed by enforcing all the larger scales in Fourier space. In contrast, in our configuration, the cloaked layer removes all the streamwise and spanwise scales, in addition to the thickness 
$l_{y,c}^+ \approx 140$ in the wall-normal direction; the present results then demonstrate that all the missing scales can be synchronized from the outer observations.
3.3. Synchronization of a vertical layer
The discussion thus far has only addressed synchronization in horizontal layers, and it is expected that the orientation of the cloaked region impacts the synchronization process. Specifically, when the layer is normal to the wall and spans the height of the channel, a distinction must be made between flow-parallel and crossflow regions because of the effect of mean-flow advection in the latter configuration. In addition, the correlation length scales of the turbulent structures are different in both directions. What remains common, however, is that the cloaked regions in both scenarios include the near-wall and outer dynamics and scales of turbulence. In this section we will examine synchronization in vertical layers that are oriented along either direction, at Reynolds number 
$Re_{\tau } = 590$.
We first consider a cloaked vertical slab that spans the height of the channel, is parallel to the flow direction and has spanwise width 
$l_z$, such that 
$\varOmega _s = \{\boldsymbol x \in [0,L_x] \times [0,L_y] \times [z_0,z_0 + l_z] \}$. The temporal evolution of streamwise-averaged synchronization error is shown in figure 8, when 
$l_z^+ = 38$. Although the initial error is uniformly 
$\sqrt {2}$ times the local r.m.s. fluctuations at different wall-normal locations (panel a), the error decreases more rapidly in the bulk than in the near-wall region (panel b–d). The dominance of near-wall synchronization error in panel (d) persists until synchronization is achieved. The inhomogeneous evolution of errors shown in figure 8 can be explained with reference to the significant variation in the Taylor microscales in the wall-normal direction. The spanwise Taylor microscale ranges from 
$2 (\varLambda _{z,u}^+, \varLambda _{z,v}^+, \varLambda _{z,w}^+) \approx (39, 22, 47)$ at the location of peak turbulence kinetic energy production, 
$y^+ = 12$, to 
$2 (\varLambda _{z,u}^+, \varLambda _{z,v}^+, \varLambda _{z,w}^+) \approx (132,114,155)$ at the channel centre. As discussed at the end of § 3.2, the smaller Taylor microscale near the wall is indicative of a shorter swept distance during the short lifetimes of the local Kolmogorov eddies. As a result, it is more difficult to synchronize the near-wall flow from the boundary observations. As such, the critical width for which synchronization is guaranteed is expected to be dictated by the near-wall dynamics. This view is reinforced by the identified critical width 
$l_{z,c}^+ \approx 45$, which is the order of the Taylor microscales in the near-wall region. The criterion 
$l_{z,c}^+ \lesssim 45$ ensures that the velocity observations on the spanwise boundaries can influence the entire layer during the lifetimes of the shortest surviving Kolmogorov eddies.
Figure 8. Synchronization of a vertical flow-parallel layer at 
$Re_{\tau }={590}$, when the layer width is 
$l_z^+ = 38$. The contours show the streamwise-averaged synchronization error, normalized by the local true root-mean-squared fluctuations 
$\mathcal {E}_{x}(u) / u^{\prime }_{rms}$. Results are shown for (a–d) 
$t^+ = \{0, 64, 191, 318 \}$.
Synchronization in a vertical, crossflow layer, 
$\varOmega _s = \{\boldsymbol x \in [x_0,x_0+l_x] \times [0,L_y] \times [0,L_z] \}$, is fundamentally different from the previous configurations due to the effect of mean advection. Intuition suggests that it is sufficient to prescribe a single crossflow plane of observations, akin to simulations of boundary layers with prescribed inflow and advective outflow conditions. When the inflow data are obtained from a reference simulation, the synchronization simulation will reproduce the reference trajectories to within machine precision. If the numerical algorithms for the reference and synchronization simulations differ, or if the inflow data are generated from experiments, the trajectories from the synchronization simulation are anticipated to diverge downstream at an exponential Lyapunov rate due to the chaotic nature of turbulence. While this intuition is helpful, it must be refined for the present channel-flow configuration because the system is closed, unlike parabolic boundary-layer flows.
We performed a number of synchronization simulations with different cloaked streamwise extents, 
$l_x$. Below a critical 
$l_{x,c}$ the flow synchronizes, and above it the flow trajectories diverge from the reference simulation. The most instructive case is reported in figure 9(ai–aiii), with 
$l_x = 2.7{\rm \pi}$ and the domain length is 
$L_x=3{\rm \pi}$. The errors within the cloaked region undergo cycles of amplification and decay, and the synchronization exponent undulates around zero. The first panel 9(ai) shows contours of the crossflow-averaged synchronization error 
$\mathcal {E}_{yz}(u)$ as a function of streamwise distance and time. The pattern of the contours is due to two effects. Firstly, the initial errors which are uniform within the cloaked region amplify exponentially as the flow advects downstream (slanted lines in the 
$(x,t)$ plane, emanating at 
$t=0$). Secondly, an influx of accurate velocity data (blue region near 
$x=0$) reduces the magnitude of the errors within the volume. The outcome is an alternation of high- and then low-magnitude errors near 
$x=l_x$, which result in an alternating pattern of high and low inflow errors at 
$x=0$, and the process repeats. The exponential amplification of errors is shown in panel (aiii), along lines of constant speed 
$\mathcal {U}=0.61$. The impact of the errors at 
$x=l_x$ on 
$x=0$ is despite the direct substitution 
$\boldsymbol {u}_s=\boldsymbol {u}_r$ within 
$\varOmega _f$, because the pressure is not observed. The pressure errors are reported in panel (aii) vs 
$x\in [0,L_x]$, and indeed its elliptic nature is evident communicating the errors across the entire length of the channel. These results suggest that a larger simulation domain, specifically 
$\varOmega _f$, could promote synchronization for this value of 
$l_x=2.7{\rm \pi}$. Panels (bi–bii) report on extending the channel from 
$L_x=3{\rm \pi}$ to 
$L_x=4{\rm \pi}$, and repeating the synchronization experiment. Indeed the cloaked layer now becomes stable and synchronizes to the reference system, both in terms of the velocity (panel bi) and pressure (panel bii).
Figure 9. Synchronization of a vertical, crossflow layer at 
$Re_{\tau }={590}$. The streamwise extent of the layer is 
$l_x = 2.7{\rm \pi}$, and two channel lengths are considered: (a) 
$L_x = 3{\rm \pi}$; (b) 
$L_x = 4{\rm \pi}$. (ai, bi) Space–time evolution of the synchronization errors in the streamwise velocity, averaged in 
$(y, z)$ and normalized by the initial value, 
$\log _{10}(\mathcal {E}_{yz}(u) / \mathcal {E}_{yz,0}(u))$. (aii, bii) Space–time evolution of error in pressure 
$\log _{10}(\mathcal {E}_{yz}(p) / \mathcal {E}_{yz,0}(u))$. (aiii) Time dependence (black) of the errors along the lines 
$x = x_0 + \mathcal {U} (t - t_0)$ and Lyapunov amplification (red dashed).
For a closed system, when the domain size is fixed, for example, circular Couette flow, increasing 
$l_{x}$ simultaneously reduces the forcing region 
$\varOmega _f$ and, hence, the critical 
$l_{x,c}$ is unambiguous. For the channel configuration, the domain length can impact the critical length for synchronization 
$l_{x,c}$. Nonetheless, the results showed that 
$l_{x,c}$ is much larger than the Taylor microscale, and instead scales with the distance travelled by the inflow during Lyapunov time scale, 
$l_{x,c} \sim O(\mathcal {U} \tau _{\sigma })$, where 
$\mathcal {U}$ is the advection velocity, 
$\tau _{\sigma } = 1 / \sigma$ is the Lyapunov time scale and 
$\sigma$ is the leading Lyapunov exponent. This criterion can also be interpreted in terms of the inverse cascade of errors (Leith & Kraichnan Reference Leith and Kraichnan1972; Boffetta & Musacchio Reference Boffetta and Musacchio2001): within 
$\varOmega _s$, perturbations at the Kolmogorov scale 
$\eta$ will contaminate the dynamics at the larger scale 
$\ell \sim 2\eta$, and the expected time of this process is proportional to the Kolmogorov time scale 
$\tau _{\eta }$; likewise, errors at 
$\ell \sim 2\eta$ will be transported towards the next larger scale 
$\ell \sim 4\eta$, and so forth. The summation of eddy turnover times at all scales is finite, and approximately equal to the Lyapunov time scale 
$\tau _{\sigma }$ (Berera & Ho Reference Berera and Ho2018). In other words, the criterion 
$l_x \lesssim O(\mathcal {U} \tau _{\sigma })$ prevents the small-scale errors from contaminating the large scales before the flow reaches the downstream boundary of 
$\varOmega _s$. The impact of the periodic boundary condition can be lessened by extending the forcing region 
$\varOmega _f$, and in the limit of very long 
$L_x$ the errors at 
$x=l_x$ inappreciably influence the inflow at 
$x=0$. The set-up then resembles predicting downstream evolution from accurate upstream data, and any infinitesimal errors due to differences in the numerical set-up still amplify due to chaos, and, hence, trajectories will diverge relative to the true system at the Lyapunov rate.
3.4. Synchronization of multiple layers
In the previous sections only one layer was cloaked and the rest of the domain was observed in order to achieve synchronization. An intuitive inquiry is the feasibility of synchronization in two or multiple layers that are individually below their respective critical thicknesses, and that are separated from each other by some observations. In this section we will examine synchronization in two horizontal layers, and then proceed to explore synchronization in multiple vertical, flow-parallel layers spaced along the spanwise directions.
A schematic for simultaneously removing two horizontal layers 
$\varOmega _{s,1} = \{\boldsymbol {x}^+ \in [0,L_x^+] \times [0,28] \times [0,L_z^+]\}$ and 
$\varOmega _{s,2} = \{\boldsymbol {x}^+ \in [0,L_x^+] \times [29,57] \times [0,L_z^+]\}$ is shown in figure 10(ai). We recall that these two layers, if removed independently, have different synchronization exponents that we can denote 
$\alpha _1$ and 
$\alpha _2$. The synchronization errors within 
$\varOmega _{s,1}$ and 
$\varOmega _{s,2}$ are plotted in panel 10(aii) when both layers are simultaneously cloaked (solid lines) and are compared with the previous experiments where only one of the layers was removed (dashed lines). The first point to note is that synchronization is successful, even though only one plane of velocity data is prescribed between the two regions. Secondly, beyond an initial transient, the error-decay rates in both layers are similar, and lie between 
$\alpha _1$ and 
$\alpha _2$. This behaviour of the errors demonstrates the co-dependence of the two layers: they do not synchronize independently of one another. Equation (3.4) governs the synchronization error averaged over the two layers, and we can derive the following expression, similar to (3.5), for their shared synchronization rate:
$$\begin{align} \alpha &= \frac{1}{\langle \| \boldsymbol e\|^2\rangle_{\varOmega_s}} \sum_{i=1,2} \frac{\varOmega_{s,i}}{\varOmega_s}(-\langle\boldsymbol{e}, \boldsymbol{e} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{r}\rangle_{\varOmega_{s,i}} - \nu \langle \| \boldsymbol{\nabla} \boldsymbol e\|^2\rangle_{\varOmega_{s,i}} + B_i )\nonumber\\ &\equiv \sum_{i=1,2} \frac{\varOmega_{s,i}}{\varOmega_s} \left(\mathcal{P}_{s,i} - \mathcal{D}_{s,i} + \mathcal{B}_{s,i} \right). \end{align}$$
This rate is determined by the net sum of four terms on the right-hand side of (3.9): production 
$\mathcal {P}_{s,i}$ and dissipation 
$\mathcal {D}_{s,i}$ in layers 
$i=\{1,2\}$. These terms depend on local quantities and therefore their magnitudes do not change appreciably whether one of two layers are cloaked. The net outcome is therefore expected to be between the two rates from the separate layers. Specifically, the more negative balance (
$\mathcal {P}_{s,2} - \mathcal {D}_{s,2}$) of the wall-detached layer promotes faster decay of errors for the entire system than when the wall layer alone is synchronized. Therefore, the asymptotic decay rate 
$\alpha$ when both layers are removed is bounded by 
$\alpha _1$ and 
$\alpha _2$
Figure 10. Synchronization of (a) two simultaneously cloaked horizontal layers and (b) multiple vertical flow-parallel layers, at 
$Re_{\tau }=590$. (ai) Schematic of the two unobserved horizontal layers. The lower layer is attached at the wall, and separated from the upper layer by one observation plane. (aii) Temporal evolution of the volume-averaged synchronization errors. Dashed lines, only one horizontal layer is cloaked: (black) 
$\varOmega _s = \varOmega _{s,1}$; (green) 
$\varOmega _s = \varOmega _{s,2}$. Solid lines: both layers are simultaneously cloaked, and errors are averaged within (black) 
$\varOmega _{s,1}$ and (green) 
$\varOmega _{s,2}$. (bi) Schematic of the cloaked, vertical, flow-parallel layers separated by one observation plane. (bii) Temporal evolution of the volume-averaged synchronization errors for 
$l_{z}^+ = 29$ (dashed line), 
$l_{z}^+ = 39$ (dashed-dotted line), and 
$l_{z}^+ = 48$ (solid line).
For synchronization in multiple layers, we examine cloaking of vertical slabs that span the height of the channel and are parallel to the flow direction, separated by observation data in the span (figure 10bi). The thickness of each unknown layer is 
$l_z$, and only one plane of observation data 
$\boldsymbol {u}_r$ is enforced between adjacent layers. Although we set the starting location of the first layer at 
$z=0$, this parameter does not affect synchronization due to the statistical homogeneity of turbulence in the span. The computational set-up is slightly different from previous sections: for both the reference and synchronization simulations, the spanwise resolution is refined to 
$\Delta z^+ = 2.4$ to ensure that each layer is well resolved. In addition, the direct substitution step 
$\boldsymbol {u}_s(t+\Delta t)=\boldsymbol {u}_r(t+\Delta t)$ in algorithm 1 is replaced by a relaxation equation 
$\boldsymbol u_s(t+\Delta t) \leftarrow (1 - \gamma )\boldsymbol u_r(t+\Delta t) + \gamma \boldsymbol u_s(t+\Delta t)$ with 
$\gamma = 0.5$, in order to ensure stability of the synchronization simulation. Additional numerical experiments using different relaxation factors, down to 
$\gamma = 0.2$, were performed and the synchronization exponents remain essentially unchanged.
The evolution of synchronization errors is reported in figure 10(bii) for 
$l_z^+ = \{29, 39, 48\}$. Similar to synchronization in horizontal layers, the errors in figure 10(bii) decay more slowly as the removed layers become thicker. The critical width below which synchronization in all the layers is guaranteed is 
$l_{z,c}^+ \approx 42$, which is slightly smaller than the criterion for removing a single layer, 
$l_{z,c}^+ \approx 45$ (cf. § 3.3). Since the balance of production and dissipation of errors per unit width is statistically stationary in the span, the difference in the critical width is due to the increase in the number of interfaces. Specifically, the cumulative contribution from the boundary terms in (3.9) increases as the number of interfaces between 
$\varOmega _s$ and 
$\varOmega _f$ increases.
The synchronization process when 
$l_z^+ = 39$ is examined in spectral space in figure 11. Since the equi-spaced observation planes can capture scales larger than 
$2 l_z^+$, the spectra at 
$\lambda _z^+ > 2 l_z^+$ are already accurately estimated near the initial time (panels (i) and green curves in panels b,d), including all the energy-containing scales and part of the dissipative scales. By contrast, the energy and enstrophy in the smaller scales 
$\lambda _z < 2 l_z$ are poorly estimated and oscillate due to the discontinuities across the observation planes. At later times (panels (ii) and blue curves in panels b,d), a wider range of scales becomes better reconstructed. Eventually the estimated spectra (colours in panels iii) converge to the true state (lines in panels iii) at all the scales, and synchronization is successfully achieved.
Figure 11. Spanwise pre-multiplied spectra of (a,b) the streamwise velocity 
$\log _{10} \langle k_z^+ |\hat {u}^+|^2 \rangle _x$ and (c,d) enstrophy 
$\log _{10} \langle k_z^+ \|\hat {\boldsymbol {\omega }}^+ \|^2 \rangle _x$, for synchronization in multiple spanwise layers with 
$l_z^+ = 39$, at 
$Re_{\tau } = 590$. The spectra are averaged in the streamwise direction, and reported at 
$t^+ = \{13,318,1270\}$. Colour contours are for (a,c) the synchronization field, and line contours at 
$t^+=1270$ are from the reference simulations. (b,d) The pre-multiplied spectra from the synchronization simulation, extracted at 
$y=1$. Black dashed lines: 
$\lambda _z^+ = 2 l_z^+ = 78$.
This configuration can also be discussed in relation to earlier studies of synchronization in spectral space. The availability of data at constant spanwise intervals is comparable to observing spanwise wavenumbers 
$k_z^+ \le 2 {\rm \pi}/ 2 l_z^+$. Successful synchronization below 
$l_{z,c}^+=42$, or equivalently above 
$k_{z,c}^+=0.0748$, is not straightforward to express in terms of the Kolmogorov length scale which depends on the wall-normal distance. For the benefit of this discussion, we remark the Kolmogorov scale is approximately 
$\eta ^+=1.66$ at the location of peak turbulence-kinetic-energy production, where the Taylor microscale was commensurate with 
$l_{z,c}^+$. Therefore, 
$l_{z,c}^+\approx 25 \eta ^+$, or 
$k_{z,c} \eta \approx 0.12$, is similar to the criterion for synchronization in spectral space (Yoshida et al. Reference Yoshida, Yamaguchi and Kaneda2005; Lalescu et al. Reference Lalescu, Meneveau and Eyink2013). Despite this favourable comparison, we note that the present configuration is different, not least because all the streamwise and wall-normal wavenumbers are unknown within the cloaked region. More importantly, within the cloaked region, the present configuration features dynamics that are unique to wall-bounded turbulence that must be synchronized within the cloaked region (e.g. wall-vorticity flux, wall-normal separation of dissipation and production, wall-normal dependence of ejections and sweeps, 
$\ldots$ , etc).
The success of synchronization of multiple layers, with merely planar observations in between layers, raises an interesting question: Given planar observations, is it possible to accurately reconstruct the velocity field in a subdomain bounded by that surface without simulating the entire system? This question is addressed in the next section.
3.5. Synchronization in subdomain simulations
A large volume of velocity observations could be difficult to acquire and, in practice, is unlikely to span the entire system beyond 
$\varOmega _s$. Even if observations are available for the complete state outside 
$\varOmega _s$, the computational cost of synchronization simulations in the full system can be prohibitive at high Reynolds numbers. It is therefore desirable to attempt synchronization of the cloaked region by simulating a subdomain of the full system and using planar time-dependent observations, which in the context of detailed laboratory experiments can be obtained using time-resolved particle image velocimetry (PIV) (Hong et al. Reference Hong, Katz, Meneveau and Schultz2012). An example using Kolmogorov flow is provided in Appendix B, where we perform synchronization simulations in a subvolume using velocity observations on the bounding planes only. In this section we retain our focus on channel flow, and return to the case of synchronization of a wall-attached layer but with an observer system that is a subvolume of the reference configuration and more limited observations.
The reference system 
$\boldsymbol {u}_r$ is still turbulent flow in the entire channel, and the velocity field is observed at one horizontal plane only, 
$\varOmega _f = \{ y = \tilde l_y \}$. The numerical scheme of the synchronization simulation 
$\tilde {\boldsymbol u}_s$ is slightly different from the reference simulation. The full Navier–Stokes equations are only solved within the wall-attached layer 
$\varOmega = \{y \in [0,\tilde l_y] \}$ instead of the entire channel. In addition, the observed velocity is enforced at 
$y=\tilde {l}_y$ as the top boundary condition, along with a homogeneous Neumann condition 
$\partial p / \partial n = 0$ on the pressure because 
$p$ is not observed. Together with periodic boundary conditions in the horizontal directions and no-slip at the wall, these six faces bound 
$\varOmega _s$, without any additional data. The initial estimate of 
$\tilde {\boldsymbol {u}}_s$ is either the mean flow superposed with white noise or the slightly disturbed true state, as explained in § 2. Our objective is to synchronize the flow in 
$\varOmega _s = \{\boldsymbol x \in [0,L_x] \times [0,\tilde l_y) \times [0,L_z]\}$ to the reference state 
$\boldsymbol u_r$.
The subdomain synchronization process at 
$Re_{\tau }=590$ with 
$\tilde l_y^+ = 28$ is visualized in figure 12(a), and is compared with the true state in panels (b). The white noise at the initial time (panel ai) decays quickly and the vortical structures near the top boundary (panel aii) become similar to the true field, although the wall-attached vortices are absent. After a sufficiently long time, all of 
$\varOmega _s$ is affected by the top boundary condition and the estimated state within the subdomain simulation synchronizes to the reference state: all the scales and structures in panel (aiii) identically match those in (
$b$iii), to within machine precision. These qualitative properties are similar to the full-domain synchronization results in figure 2, although here the observer system is only a subdomain simulation with height 
$l_y^+=28$.
Figure 12. Instantaneous vortical structures visualized using the 
$\lambda _2$ vortex identification criterion with threshold 
$\lambda _2 = -4$. (a) Synchronization simulation in a wall layer 
$l_y^+ = 28$. (b) Reference simulation. Results are shown for (i–iii) 
$t^+= \{0, 13, 1000\}$.
The errors from the subdomain synchronization are reported in figure 13(a) (blue lines), where they are compared with results using the approach from § 3.1 for full-domain synchronization (black lines). We verified that the critical layer thickness remains unchanged, 
$\tilde {l}^+_{y,c} = l^+_{y,c} \approx 32$. When the subdomain height is larger than this value, 
$\tilde {l}_y > l_{y,c}$, the errors diverge similar to the previously reported behaviour in § 3.1. For smaller layer heights, below the critical value, the errors decay exponentially which indicates convergence towards the true flow trajectories, initially at the same rate as the full-system configuration. However, a significant difference can be observed: the errors in the present subdomain approach do not decay to machine precision, but rather saturate at approximately 
$10^{-3}$ of the initial value, or equivalently, three orders of magnitude lower than the reference r.m.s. fluctuations. These persistent long-time errors are not due to failure of synchronization of the dynamics as described by the Navier–Stokes equations, but are rather due to the differences in the numerical solutions of the governing equations in a subdomain vs the full channel, specifically the fractional-step method. In support of this argument, we performed a new reference simulation 
$\tilde {\boldsymbol u}_r$ within the wall-attached layer, i.e. same subdomain as the synchronization simulation and with the same numerical method. The only difference between 
$\tilde {\boldsymbol {u}}_s$ and 
$\tilde {\boldsymbol {u}}_r$ is the initial condition, the former starting from white noise and the latter from the true initial condition. This new reference simulation 
$\tilde {\boldsymbol u}_r$ does not match the original 
$\boldsymbol {u}_r$ from the global domain, and the difference is of the order of 
$10^{-3}$ – a matter that we return to below. More importantly, the synchronization solution 
$\tilde {\boldsymbol u}_s$ from an arbitrary initial condition converges to 
$\tilde {\boldsymbol u}_r$ to within machine precision (figure 13b); the subsystem is therefore Lyapunov stable and can accurately synchronize when 
$\tilde {l}_{y} < \tilde {l}_{y,c} = l_{y,c}$.
Figure 13. Synchronization in a subdomain of the channel, at 
$Re_{\tau }=590$. (a) Volume-averaged errors 
$\mathcal {E}_{s}$ of the instantaneous velocity, when synchronization simulations are performed in a subdomain (blue), compared with the full domain (black). The reference velocity 
$\boldsymbol u_r$ in the definition of the errors was obtained from a full-domain simulation. (b) The errors of the subdomain synchronization experiments are re-evaluated using a reference velocity from a simulation in the same subdomain, 
$\tilde {\mathcal {E}}_{s} = \langle \|\tilde {\boldsymbol u}_s - \tilde {\boldsymbol u}_r \|^2\rangle _{\varOmega _s}^{1/2}$. Dashed lines, 
$\tilde l_y^+ =l_y^+ = 18$; dashed-dotted lines, 
$\tilde l_y^+= l_y^+=28$; solid lines, 
$\tilde l_y^+= l_y^+=38$.
The reason that the reference simulations in the sub- and full domains, 
$\tilde {\boldsymbol u}_r$ and 
$\boldsymbol u_r$, differ can be explained using a similar analysis as Wu, Zaki & Meneveau (Reference Wu, Zaki and Meneveau2020). Briefly, in addition to the same initial condition, the intermediate velocity of the fractional-step algorithm for solving the Navier–Stokes equations must be enforced at 
$y = \tilde {l}_y$ to eliminate the mismatch. This velocity is not physical, and is only introduced for the purpose of numerically solving the equations. We have demonstrated synchronization to 
$\tilde {\boldsymbol u}_r$ in the subsystem, and the deviation from the full system is a matter of a difference in the numerical model. Notwithstanding this technical detail, the results presented here and in Appendix B demonstrate that with only planar observations, synchronization in a subdomain that satisfies the critical size criteria can accurately converge onto the true trajectories of the full system, to within the errors of the numerical model (figure 13a,b), and, hence, is sufficiently accurate to evaluate any flow quantity of interest including velocity gradients or vortical structures (cf. figure 12).
4. Conclusions
Synchronization of chaos in turbulent channel flow was attempted using continuous data assimilation techniques. Fully resolved observations of the turbulence are available outside a cloaked, or unobserved, region of the flow. These observations are directly substituted in the synchronization simulation at every time step in order to drive the missing flow towards the reference state. The temporal evolution of the estimation error was adopted to evaluate a synchronization exponent, which determines the success and rate of synchronization. When synchronization is successful, all the turbulence scales in the cloaked region are accurately re-established, and the synchronization exponent is independent of the initial condition.
Synchronization was examined in detail for a cloaked, wall-attached horizontal layer. After the initial transient, a successful synchronization features monotonic decay of the errors within the layer. The synchronization rate, when scaled in wall units, is independent of the instantaneous reference state and the Reynolds number up to 
$Re_{\tau } = 1000$. The critical layer thickness is 
$l_{y,c}^+ \approx 32$, which indicates that the near-wall turbulence, including the true vorticity source, are interpretable from outer observations.
When the unknown horizontal layer is detached from the wall, its synchronization rate and critical thickness depend on its wall-normal distance, 
$y_0$. A thin layer has a relatively slow synchronization rate compared with a wall-attached counterpart, but the critical thickness increases monotonically with 
$y_0$. Up to 
$y_0^+ = 100$, the synchronization rate and critical thickness expressed in wall units are independent of the Reynolds number. These trends are the outcome of the balance in the equation governing the synchronization errors, which features a source term at the boundary of the cloaked region, production of errors against the mean shear and viscous dissipation of errors. Relative to the scales of turbulence, the critical thickness of a cloaked layer is comparable to twice the Taylor microscale which quantifies the distance traversed by Kolmogorov eddies within their lifetimes as they are transported by the r.m.s. velocity.
Simultaneous synchronization in multiple volumes was investigated by considering two horizontal layers separated by a single plane of observations. The synchronization rate was bounded by the exponents that correspond to cloaking each layer alone. We also examined synchronization of vertical layers that are parallel to the flow, and separated by a single observation plane. The critical layer width that ensured synchronization was identified, 
$l_{z,c}^+ \approx 42$, and is of the order of the near-wall Taylor microscale of the turbulence. The criterion for synchronization of a vertical, crossflow volume is set by the mean advection and the Lyapunov time scale, and is therefore less restrictive than the criteria in the other directions for the range of Reynolds numbers examined herein.
The present results thus extend the synchronization criteria for isotropic turbulence to wall-bounded turbulence. The identified conditions based on the Taylor microscale and Lyapunov time scale accommodate inhomogeneity and anisotropy of turbulence and, therefore, can be applied in various configurations.
Our final experiments, which are of practical interest, attempted reconstruction of the cloaked region using simulations of a subdomain of the reference system, and with limited observations. We revisited synchronization of a wall layer, when the available data are one plane of instantaneous velocity observed from the reference channel-flow system. The observer system is then the cloaked near-wall layer and that plane of data only. The reconstructed flow fields perfectly matched resimulations of the subdomain using the true initial conditions, and both exhibited some deviation from the true, full channel-flow system due to the discretization scheme of the Navier–Stokes equations. This mismatch can be reduced by increasing the number of observed velocity planes. The practical implication of these results is the possibility of accurately predicting the near-wall flow and the wall stresses from velocity measurements at a wall-parallel plane at 
$y^+ \lesssim 32$; if the measurements are performed beyond 
$y^+ = 32$, variational data assimilation techniques must be adopted and can provide bounds on the accuracy of the estimated near-wall flow as a function of the measurement height. A more general configuration was investigated in Appendix B for Kolmogorov flow. We demonstrated the feasibility of synchronization in a subdomain using observations on the surrounding six faces only, as long as the critical dimensions are obeyed. All together, the present results demonstrate the potential of augmenting detailed experimental observations, e.g. time-resolved planar PIV measurements, using continuous data assimilation techniques that are efficient, easy to implement, and can provide access to entire true flow trajectories.
Funding
This work was supported in part by the Office of Naval Research (grant N00014-20-1-2715). Computational resources were provided by the Maryland Advanced Research Computing Center (MARCC).
Declaration of interests
The authors report no conflict of interest.
Appendix A. Equation for volume-averaged synchronization error
In this appendix we outline the procedures to derive the equation for the volume-averaged synchronization error, 
$\mathcal {E}_{s} = \langle \|\boldsymbol e \|^2 \rangle _{\varOmega _s}^{1/2}$. By subtracting the Navier–Stokes equations for 
$\boldsymbol u_s$ and 
$\boldsymbol u_r$, we obtain the equations governing 
$\boldsymbol e := \boldsymbol u_s - \boldsymbol u_r$ and 
$\zeta := p_s - p_r$,
\begin{gather} \frac{\partial \boldsymbol{e}}{\partial t}+\boldsymbol{e} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{r}+\boldsymbol{u}_{r} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{e}+\boldsymbol{e} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{e} ={-}\boldsymbol{\nabla} \zeta+\nu \nabla^{2} \boldsymbol{e} + \boldsymbol f, \end{gather}
\begin{gather}\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol e = b(\boldsymbol x_b) \delta(\boldsymbol{x} - \boldsymbol{x}_{b}). \end{gather}
The forcing 
$\boldsymbol f$ represents direct substitution of the observations in 
$\varOmega _f$, which ensures that 
$\boldsymbol e = \boldsymbol 0$ for all 
$\boldsymbol {x} \in \varOmega _f$. Note that the pressure difference 
$\zeta$ is non-zero through the entire channel because the reference pressure 
$p_r$ is not observed. Due to discontinuity between the observed and estimated velocity fields at all 
$\boldsymbol x_b \in \partial \varOmega _s$, the corresponding divergence of the synchronization error is not zero, which leads to the source term 
$b(\boldsymbol x_b)$ in (A1). The specific form of the source term depends on the numerical discretization. By calculating the dot product of (A1) with 
$\boldsymbol e$ and then averaging over 
$\varOmega _s$, the equation for 
$\mathcal {E}_{s}^2 = \langle \|\boldsymbol e \|^2 \rangle _{\varOmega _s}$ can be derived, i.e.
\begin{equation} \frac{1}{2} \frac{{\rm d}}{{\rm d} t}\langle \| \boldsymbol e\|^2\rangle_{\varOmega_s} ={-}\langle\boldsymbol{e}, \boldsymbol{e} \boldsymbol{\cdot} \boldsymbol{\nabla} \boldsymbol{u}_{r}\rangle_{\varOmega_s} -\nu \langle \| \boldsymbol{\nabla} \boldsymbol e\|^2\rangle_{\varOmega_s}+B. \end{equation}
Other terms that can be written in conservative form vanish, either by virtue of 
$\boldsymbol e = \boldsymbol 0$ on the boundaries of 
$\varOmega _s$ or appropriate boundary conditions, for example, periodicity in the horizontal directions for a horizontal slab and no-slip conditions on the wall. The boundary term in the kinetic energy equation is
\begin{equation} B = \frac{1}{\varOmega_s} \int_{\varOmega_s} b(\boldsymbol x_b) \delta (\boldsymbol x - \boldsymbol x_b) \zeta(\boldsymbol x)\,{\rm d}V = \frac{1}{\varOmega_s} \int_{\partial \varOmega_s} b(\boldsymbol x_b) \zeta(\boldsymbol x_b) \,{\rm d}S. \end{equation}
This term only contains a surface integration over the boundary 
$\partial \varOmega _s$ and is generally much smaller than the production and dissipation in (A3). For the wall-attached horizontal layer 
$\varOmega _s = \{\boldsymbol x \in [0,L_x] \times [0,l_y] \times [0,L_z]\}$ discussed in § 3.1, the boundary term (A4) becomes
\begin{equation} B = \frac{1}{\varOmega_s} \int_{\varOmega_s} b(x,z) \delta (y-l_y) \zeta(x,y,z)\,{\rm d}V = \frac{1}{\varOmega_s} \int_{\partial \varOmega_s} b(x,z) \zeta(x,l_y,z) \,{\rm d}\kern0.06em x\,{\rm d}z. \end{equation}Appendix B. Subdomain synchronization in Kolmogorov flow
In § 3.5 we demonstrated that the wall layer in turbulent channel flow can be accurately recovered by performing the synchronization simulation in a subdomain, with the limited observations of the velocity on the top boundary only. The synchronization subdomain thus spanned the horizontal dimensions of the entire channel, with all the streamwise and spanwise turbulent scales being unknown. Motivated by practical limitations, in this appendix we consider the more specialized case where the synchronization subdomain is a general, three-dimensional, rectangular subvolume that is neither periodic nor attached to a boundary, and the observations are limited to the velocity data on the surrounding planes only. The synchronization experiments are performed in Kolmogorov flow which was the subject of previous studies of synchronization, albeit solely in Fourier space (Lalescu et al. Reference Lalescu, Meneveau and Eyink2013).
The reference simulation is performed in a tri-periodic, cubic box with sides 
$L = 2{\rm \pi}$, and driven by a body force 
$\boldsymbol f_b = (0.2 \sin y, 0, 0)$. The Reynolds number is 
$Re_{\varLambda } =u^{\prime }_{rms} \varLambda /\nu = 146$, where the Taylor microscale 
$\varLambda = \sqrt {15 \nu / \mathcal {D}}\, u^{\prime }_{rms}$ is evaluated using global dissipation 
$\mathcal {D}$. The computational domain is discretized on a uniform grid with 256 points in each direction. The incompressible Navier–Stokes equations are solved using the same numerical algorithm as in § 2. Due to periodicity, the pressure Poisson equation is solved using Fourier transform in all three directions. Once the flow reaches a statistically stationary state (see sample snapshots in figure 14a), instantaneous velocity data on six faces of a cubic subdomain 
$[l_0, l_0 + \tilde l]^3$ are extracted and stored.
Figure 14. Instantaneous visualizations of synchronization in Kolmogorov flow. (a) Contours of the reference spanwise velocity 
$w_r$, plotted at 
$x=L/2$. (b) Enlarged, three-dimensional views of the reference subdomain of interest and (c) of synchronization using a subdomain simulation. (a–c) The three snapshots from bottom to top correspond to 
$t=\{0, 1, 4\}$.
The synchronization simulations are performed in subdomains bounded by these six faces, and the velocity observations are prescribed as Dirichlet boundary conditions. Since pressure is not observed, Neumann boundary conditions 
$\partial p / \partial n = 0$ are adopted, and Fourier-cosine transforms are performed to solve the pressure Poisson equation. The starting location of the subdomain is fixed 
$l_0 = L/16$, and three subdomain sizes 
$\tilde {l}=\{{1}/{8}, {1}/{4}, {1}/{2}\}L$ are considered. The initial estimate for the first two cases is white noise proportional to the root-mean-squared velocity fluctuations, and a Lyapunov-type experiment is performed for 
$\tilde l=L/2$.
The subdomain synchronization simulation for 
$\tilde l = L/4$ is visualized in figure 14(c) at 
$t=\{0, 1, 4\}$, and compared with the true state in panel (b). The velocity field on the surface of the cube is identical to the true state because the continuous data assimilation in this configuration is effectively prescribed as the time-dependent velocity boundary condition. Deviations from the reference state occur only at interior points (see figure 14b,c), and diminish as the flow evolves in time. The volume-averaged synchronization error is plotted in figure 15 (dashed-dotted line), normalized by its initial value. The error decays to 
$10^{-4}$ and then saturates, similar to the behaviour discussed in § 3.5 when the reference velocity is obtained from the full-domain simulation. Note that the error would decay to machine precision if compared with a reference simulation in the subdomain, starting from an accurate initial condition.
Figure 15. Temporal behaviour of errors in subdomain synchronization simulations of Kolmogorov flow. The volume-averaged errors 
$\mathcal {E}_{s}$ are normalized by their initial value 
$\mathcal {E}_{s,0}$. Dashed line, 
$\tilde l = L/8$; dashed-dotted line, 
$\tilde l = L/4$; solid line, 
$\tilde l = L/2$.
Figure 15 shows the synchronization errors for all three cases, 
$\tilde l=\{{1}/{8}, {1}/{4}, {1}/{2}\}L$. The critical subdomain size 
$\tilde {l}_c$ lies between 
$\tilde {l}=L/4$ and 
$\tilde {l}=L/2$. Because the large-scale flow is inhomogeneous in the 
$y$ direction, the local Taylor microscales are evaluated at each 
$y$ position from the two-point correlations. The maximum value over 
$y$ and the 
$\{u ,v, w\}$ components are 
$2\varLambda _{x,u}=0.95$, 
$2\varLambda _{y,v}=0.82$ and 
$2\varLambda _{z,w}=0.75$, which are all smaller than 
$L/4$. Similar to removing streamwise layers in § 3.4, the critical subdomain size in 
$x$ exceeds 
$2\varLambda _x$ because of the mean advection effect. In effect, as long as the subdomain size is below the threshold for synchronization in one direction, the flow can synchronize to the reference state by enforcing planar velocity data on the six faces.
Open access
