{'fraction_non_alphanumeric': 0.03007642100429663, 'fraction_numerical': 0.008101282429728754, 'mean_word_length': 4.847621787495205, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Model-based multi-sensor fusion for reconstructing wall-bounded turbulence

Wall-bounded turbulent flows can be challenging to measure within experiments due to the breadth of spatial and temporal scales inherent in such flows. Instrumentation capable of obtaining time-resolved data (e.g., Hot-Wire Anemometers) tends to be restricted to spatially-localized point measurements; likewise, instrumentation capable of achieving spatially-resolved field measurements (e.g., Particle Image Velocimetry) tends to lack the sampling rates needed to attain time-resolution in many such flows. In this study, we propose to fuse measurements from multi-rate and multi-fidelity sensors with predictions from a physics-based model to reconstruct the spatiotemporal evolution of a wall-bounded turbulent flow. A"fast"filter is formulated to assimilate high-rate point measurements with estimates from a linear model derived from the Navier-Stokes equations. Additionally, a"slow"filter is used to update the reconstruction every time a new field measurement becomes available. By marching through the data both forward and backward in time, we are able to reconstruct the turbulent flow with greater spatiotemporal resolution than either sensing modality alone. We demonstrate the approach using direct numerical simulations of a turbulent channel flow from the Johns Hopkins Turbulence Database. A statistical analysis of the model-based multi-sensor fusion approach is also conducted.


Introduction
An ability to predict and control turbulence would benefit numerous engineering and scientific applications. Yet, such efforts are often compromised by the sensing and diagnostic requirements of resolving the breadth of spatiotemporal scales underlying turbulence. Instruments capable of resolving the fastest velocity fluctuations in time-e.g., Hot-Wire Anemometers (HWA), Laser-Doppler Anemometers (LDA)are often restricted to localized measurements at a single point in space. On the other hand, Particle Image Velocimetry (PIV) systems capable of acquiring spatiallyresolved field measurements tend to lack the temporal resolution required for higher Reynolds number flows. In this study, we formulate a framework for reconstructing the spatiotemporal evolution of turbulent flows by fusing multi-rate and multi-fidelity measurements from these different sensing modalities together with predictions from simple physics-based models.
The notion of flow estimation and reconstruction from available measurements is not a new one. Numerous efforts have considered linear stochastic estimation (LSE) to extract coherent structures from various measurement sources [1,2,10]. For example, LSE has been used in the investigation of wall-bounded turbulence and cavity flow using wall-based pressure and shear stress measurements [6,21,22]. A Proper Orthogonal Decomposition (POD)-based LSE has been utilized to address challenges of multi-scale, nonlinear, and higher order estimation [4,20,23]. Time-delay variants of the LSE-POD approaches have been investigated for real-time flow estimation [5] and flow control [26,29]. Recently, LSE and spectral POD algorithms have been combined into a unified framework that enables modal analysis of fluid flows from nontime-resolved PIV data in conjunction with unsteady pressure measurements [32]. Similarly, modal analysis based on the resolvent framework has been shown to be suitable for estimating space-time statistics of flows [27] and also useful for estimation of turbulent flows from wall-measurements [3].
Some recent investigations have considered machine learning methods and superresolution data reconstruction techniques for reproducing turbulent flows [7,8,13]. Once trained on high-quality data, these methods have been shown capable of reproducing the underlying flow field using remarkably coarse measurements. Previous studies have also considered purely statistical data fusion methods for flow reconstruction. In [30], a "model-free" maximum a posteriori (MAP) algorithm was proposed for fusing low-temporal-high-spatial resolution data with high-temporal-lowspatial resolution data for turbulent flow reconstruction. However, it is important to note that this work did not leverage a model to perform dynamic estimation. Indeed, model-based dynamic estimation can improve reconstruction performance because it offers an additional information source in the form of model predictions.
Dynamic estimation has been investigated in prior works on flow reconstruction, especially in the context of band-limited oscillator-type flows. Tu et al. [28] established a dynamic estimator to reconstruct the wake of a thick flat plate at low Reynolds numbers. This was done using a data-driven model of the flow in conjunction with a Kalman smoother that was used to fuse model predictions with non-timeresolved PIV and time-resolved point-sensor measurements. Recent efforts have focused on using physics-based models grounded in the resolvent formalism to predict and reconstruct flows [19,12,27]. The use of low-complexity physics-based models for flow reconstruction has the benefit of ensuring that the estimated flow satisfies physical constraints (e.g., conservation of mass and momentum). Recent studies have also shown that it is possible to perform a data-driven refinement of the linearized Navier-Stokes equations (NSE) to account for second-order statistics [14,31].
In previous work [16,17], we have shown that even simple linear time-invariant models can be used to reconstruct the spatiotemporal evolution of wall-bounded turbulence from non-time-resolved 2-dimensional 2-component (2D2C) field measurements, as would be obtained from a PIV system. For example, linear models grounded in Rapid Distortion Theory (RDT) [11,25] are capable of predicting the flow evolution over a time period much shorter than a typical eddy turnover time. This is a consequence of the fact that persistent features of the local eddy structure change slowly over the time of distortion, making higher-order nonlinear effects negligible over a fast time-scale [11]. As a result, we found that it is possible to use linear RDT models to reconstruct the time evolution of wall-bounded turbulent flows between snapshots of a non-time-resolved PIV system [16,17]. The model-based flow reconstruction was further improved by propagating the snapshots both forward and backward in time, then using an appropriately designed spatiotemporal weighting scheme to fuse the model predictions.
In this paper, we build upon our previous study on model-based turbulent flow reconstruction with an aim to improve performance even further. In particular, we investigate the potential for improving the quality of turbulent flow reconstruction by introducing additional time-resolved measurements at distinct spatial locations in the flow field. To do so, we propose a model-based multi-sensor fusion approach for wall-bounded turbulence reconstruction. The framework fuses model-predictions with time-resolved point measurements using a "fast filter" that can be applied both forward and backward in time. Predictions from the fast filter are fused with nontime-resolved field measurements on a slower time scale, corresponding to the sampling rate of the slow measurements. We show that both filtering tasks can be conducted with conventional Kalman filtering algorithms formulated in conjunction with the RDT-based dynamics. This also serves to improve the flow reconstruction based on noisy measurements. The model-based multi-sensor fusion approach is investigated using DNS data of a turbulent channel flow from the Johns Hopkins Turbulence Database (JHTDB) [9]. In order to fully evaluate the approach, we also investigate the role of (1) placement of the "fast" point sensors, and (2) design of the weighting function used for the forward-backward fusion on the resulting reconstruction.
The paper proceeds as follows: In Section 2, we present the ingredients needed for the reconstruction algorithm: the RDT-based dynamics model, the sensor measurement models, and the forward-backward filtering algorithms and corresponding weighting functions for data fusion. Section 3 describes the benchmark dataset and the metrics that will be used to evaluate performance of the reconstruction. Section 4 reports the results of Monte Carlo simulations for different fusion approaches and sensor placement strategies. Finally, we provide concluding remarks of our study in Section 5. iv Mengying Wang et al.

Model Formulation
An essential element in reconstructing the turbulent flow between two consecutive non-time-resolved snapshots is a model for approximating the evolution of the flow. Here, we use a linear model grounded in Rapid Distortion Theory (RDT), which is based on the assumption that nonlinear flow interactions can be neglected on timescales that are much shorter than a typical eddy turnover time [11]. This allows the equations of motion to be linearized for predictions of the flow evolution over short time-horizons. Beginning with the incompressible Navier-Stokes equations (NSE) and considering the dynamics of turbulent velocity fluctuations u = (u 1 , u 2 , u 3 ) about a mean flow U = (U (x 2 ), 0, 0)-see Fig. 1-we have where p represents pressure fluctuations and Re τ denotes the friction Reynolds number. By the assumptions of RDT, the nonlinear term NL can be neglected when the sampling period of the slow measurement system is less than a typical eddy turnover time [11]. For additional simplicity, we also relax the incompressibility constraint.
Since most commonly used PIV systems can only capture planar two-dimensional two-component snapshots, we further simplify this linear model by neglecting the out-of-plane flow and pressure gradient terms. Then, the simplified two-dimensional flow model will be: For the numerical simulation, Eq. 3 are numerically integrated in time using an explicit Euler method. A first-order upwinding scheme is used for the convection terms, and a second-order central differencing scheme is used for the diffusion and coupling terms. This simplified model can be used to predict the evolution of the flow both before and after a given snapshot, simply by propagating the dynamics either backward or forward in time, respectively. For instance, given a snapshot at time t = 0, numerical integration of Eq. 3 can be used to predict the forward evolution of the flow up to the time t = T when a new snapshot is available. Similarly, the flow can be propagated backward in time from t = T to t = 0 by utilizing the transformation τ = T − t in Eq. 3, such that the new time variable has value τ = 0 at the end of the interval. This transformation yields a negative convection velocity, mean shear, and effective viscosity in the backward dynamics. The negative viscous diffusion steepens spatial gradients in velocity, which eventually leads to blowup in the model prediction, especially when additional error is introduced by finite precision numerics. However, the blowup in the backward reconstruction due to this negative diffusion is not prominent for a short prediction horizon, as analyzed in a benchmark study [17]. Indeed, by propagating forward and backward, it is possible to formulate physically motivated weighting schemes to fuse forward and backward estimates of the flow field to improve accuracy in the reconstruction between PIV snapshots [17].
To facilitate model-based multi-sensor fusion, we discretize Eq. 3 in both space and time to obtain a finite-dimensional discrete-time state-space representation of the flow: where q = (u 1 , u 2 ) denotes the N -dimensional flow state, A + denotes the linear propagator for the forward dynamics, ∆t denotes the time increment, and w(t) denotes model uncertainty or process noise. We assume the process noise w(t) is zeromean, Gaussian distributed white noise with covariance Q, i.e., w ∼ N (0, Q). Although this is a crude modeling assumption, it will simplify the formulation of the filter, which can be carefully tuned to achieve adequate reconstruction performance. In addition, we define output equations to map the flow state at a given instant to the sensor measurements. Here, we are interested in modeling sensors with disparate sampling rates -a slow-in-time field measurement (e.g., PIV) with sampling period ∆t + s , and a fast-in-time point measurement (e.g., HWA) with sampling period ∆t + f , where the superscript "+" denotes the normalization with respect to friction velocity and kinematic viscosity of the turbulent flow. For simplicity, we assume the slow and fast time scales are integer multiples of each other, i.e., ∆t + s = ∆t + f where is an integer. This yields the fluid dynamics and sensor-output equations, where n is the time-index for the fast system. Additionally, y f ∈ R m denotes a vector of "fast" velocity fluctuation measurements -herein referred to as "fast" measurements, for simplicity -taken at m N discrete points in space, determined by a Fig. 2: Block diagram of the multi-output state-space system in Eq. 5 spatial sub-sampling matrix S ∈ R m×N composed of zeros and ones. These measurements are assumed to be contaminated by additive zero-mean Gaussian white noise v f ∼ N (0, R f ), where R f denotes the measurement noise covariance. Slow measurements are modeled in Eq. 5c, where y s denotes the non-time-resolved field measurements -herein referred to as "slow" measurements, for simplicity -which are also assumed to be contaminated with zero-mean Gaussian white noise with covariance R s : v s ∼ N (0, R s ). Note that "mod" denotes the remainder operator, so that the output equation Eq. 5c reflects that new field measurements are only obtained at slow sampling times. The block diagram of the discrete-time state space system is shown in Fig. 2.

Model-based multi-sensor fusion framework
With the dynamics and sensor models, we are now positioned to formulate a multisensor fusion approach. The general idea is to fuse the fast measurements with the model in Eq. 5 using a "fast filter". Then, a "slow filter" is formulated to fuse the fast filter estimates with the slow measurements. Since the model equations in Eq. 5 are linear models with zero-mean Gaussian distributed noise, a Kalman filter can be used in this fusion approach, as will be described momentarily. Further, the fast and slow filters can be applied both forward and backward in time in post-process. A schematic of the overarching fusion algorithm is shown in Fig. 3. Consider the reconstruction forward in time as an example. We use the first snapshot as the initial condition and apply the RDT model to propagate the system state. We fuse the fast measurement and the model prediction to yield an estimate from the fast filter, and when the next snapshot is available, a slow fusion algorithm is applied to integrate the fast estimate with the slow field measurement. The backward fusion method can be applied similarly with the final snapshot as the initial condition. Then, a weighted sum of the forward and backward filter estimates is used to reconstruct the flow with higher fidelity.
A filtering algorithm computes the marginally estimated distribution of the state from past measurements. The Kalman filter deals with the particular case where the dynamic model is linear with the zero-mean Gaussian distributed process and mea- shows the multi-rate measurement setup. The green meshed squares represent the slow measurements (e.g., from PIV) and the red track represents a sample signal from a fast measurement at a fixed location in the flow (e.g., from HWA). (b) and (c) show the "fast" filter within one interval forward and backward in time, respectively. (d) shows the fusion of estimates from the forward and backward filters together using a weighting function to reconstruct the flow.
surement noise [15]. The RDT-based model derived above meets these conditions, and so we will design a Kalman filter to handle multi-rate multi-fidelity measurements to improve the estimates of the flow state predicted by the RDT model. During the time period between two consecutive snapshots, we apply a fast Kalman filter to fuse the fast measurement with the model generated state. At each time a new snapshot is collected, we use a slow Kalman filter to fuse the slow measurement with the prediction of the state governed by the model. The method we use in this study can be viewed as a special case of time-varying Kalman filter, based on time-invariant dynamics and a time-varying measurement equation. We summarize the multi-rate forward Kalman filtering algorithm as follows: Slow filter, when (n mod ) = 0 Fast filter, when (n mod ) = 0 : where n denotes the new time index corresponding to the discretized version of the RDT model,q + (n|n − 1) denotes the a priori state estimate, andq + (n|n) denotes the a posteriori state estimate. The + subscripts on A + andq + indicate that these variables are related to the forward estimates. Note that both the slow Kalman filter and the fast Kalman filter estimate the same state q and its covariance matrix P , and that the dynamics of the state do not change with time, only the output equations do. Thus, the process uncertainty covariance matrix Q in both filters is consistent. However, the measurement noise covariance depends on the measurement system: Recall that R s and R f denote the measurement noise covariance matrix for slowin-time and fast-in-time measurements, respectively, and their distributions follow v f ∼ N (0, R f ) and v s ∼ N (0, R s ) in Eq. 5. For the backward Kalman filter, the a priori state estimate becomesq − (n|n + 1) = A −q− (n + 1|n + 1), i.e., the dynamic evolves backward in time, where the operator A − is the backward RDT model.
Establishing the multi-rate Kalman filter is expected to improve the flow reconstruction compared to using a filter based on a single measurement source alone. However, performance will depend on the specific weighting scheme used for the fusion of forward and backward estimates. A smoothing algorithm is usually used to reduce the large errors associated with forward and backward filtering at later times. Traditionally, a fixed-interval smoothing algorithm can be used directly. However, the process noise in our model accounts for uncertainties due to the neglected nonlinear terms, including the nonlinear interactions in NSE. As such, we lack an exact characterization of this process noise, and the conventional smoother requires careful tuning. In our previous study, a temporal weighting function and a spatiotemporal weighting function were applied to fuse the forward and backward model predictions with noteworthy performance [17]. As such, we opt to use the same approach for fusion here, as described next.
Assume the evolution of the turbulent flow can be propagated forward in time starting from the first slow-in-time snapshot as well as backward in time starting from the second snapshot. A smoothing algorithm of the forward and backward propagation can be formulated to achieve an improved estimate of the flow profile by implementing a fixed-interval smoothing algorithm. The resultant reconstruction of the flow field can be written in the form of a weighted combination of these forward and backward estimates, as follows:q whereq + andq − denotes the estimated outputs of the forward and backward Kalman filter, respectively, while G + and G − are the weighting factors to be determined. In this study, we consider two types of weights to be applied in this smoothing problem: (1) temporal weights and (2) spatiotemporal weights [17]. A temporal weighting scheme is formulated based on linear weighting in time between two consecutive slow-in-time snapshots as: Recalling that ∆t + s is the sampling time of the snapshots, which is also the interval period, we denote t ∈ [0, ∆t + s ] as the time within the interval. The temporal weights G + and G − ensure that the forward estimate is emphasized closer to the beginning snapshot and gradually yields its weight to the backward estimate towards the ending snapshot.
In addition, a spatiotemporal weighting scheme can be formulated to account for the hyperbolic nature of the governing equations. Recall that Eq. 3 contains a hyperbolic advection term, which is dominant in wall-bounded turbulent flows. The flow estimate follows a propagation at a speed corresponding to the local mean velocity U (x 2 ). This property yields two special regions in the x 1 − t plane: One is excluded from the region of influence (ROI) for the initial snapshot and the other is excluded from the domain of dependence (DOD) for the ending snapshot, as shown in Fig. 4. The characteristic lines with slope dt/dx 1 = 1/U (x 2 ) separate the ROI of snapshot 1 and the DOD of snapshot 2, and the common region in green illustrates the combination of forward and backward predictions. Note that the slope of the characteristic lines are determined by the local mean velocity U (x 2 ), which also varies with the wall-normal direction x 2 . The spatiotemporal weighting scheme accounts for these where l + = U (x 2 )t and l − = U (x 2 )(∆t + s −t) are calculated as the convective length scale corresponding to the local mean flow velocity U (x 2 ). L x1 is the streamwise length of the snapshot window. More details about these two approaches can be found in [17]. The reconstructions using these two approaches are compared in Section 4.

Data generation and evaluation metrics
In this section, we present some prerequisites to conduct the reconstruction and evaluate its performance: Section 3.1 introduces the DNS data source serving as ground truth for the algorithm validation in this study. The measurement signals in this paper are also simulated by subsampling the data and artificially adding sensor noise. Section 3.2 presents the metrics for quantifying the reconstruction performance.

DNS database
In this study, we consider the DNS data for a turbulent channel flow at Re τ = u τ h/ν = 1000 from the JHTDB [9], where u τ is the friction velocity, ν is the kinematic viscosity, and h is the half-height of the channel. We utilize the streamwise and wall-normal components of the velocity and subsample the data on a restricted domain of n x1 × n x2 = 65 × 65. The domain consists of a h × h planar window with a streamwise length h and vertically spanning from the bottom wall to the centerline of the channel. The spatial resolution for this dataset is ∆x + 1 = ∆x + 2 = h + /64 ≈ 16. This resolution is chosen to be consistent with a physical PIV system with analysis based on 16 × 16 pixel segments with 50% overlap as a benchmark case involving the use of a camera with 1 MP resolution. More details about the simulation parameters of the turbulence database can be found in [9]. The velocity data is extracted at the time step t + = 0.065, based on the DNS time step. The superscript "+" denotes the normalization with respect to u τ and ν such that the dimensionless h + = Re τ . Note that the plus mark "+" denotes normalization only when it appears as a superscript; when in the subscript, it indicates that the variables are related to the forward estimates, e.g., as in G + andq + .
To sample the slow measurements, we add zero-mean Gaussian distributed noise to every 2-dimensional 2-component snapshot on the slow-sampling interval ∆t + s = 6.24. The interval time ∆t + s is estimated for a PIV system capable of 100 Hz sampling rate (i.e, sampling time T = 0.01 s), corresponding to water flow with friction velocity u τ = ∆t + s ν/T ≈ 0.025 ms −1 . As discussed in our previous work [17], when the total time period between the first and last snapshot is ∆t + s ≥ 24, the model-based reconstructions become less accurate. We choose the time interval for reconstruction to be ∆t + s = 96t + = 6.24 to satisfy this condition. The noise level of the measurements is quantified by the Signal to Noise Ratio (SNR), which is defined as the ratio between local data magnitude to the variance of the additive noise. For slow measurements, we use SN R s = 20, which is a conservative estimate based on experimental experience [24].
The fast measurements are spatially subsampled from the DNS using the sparse sampling matrix S. The fast measurements are assumed to be of higher precision than the slow measurements, and contaminated with additive noise such that SN R f = 100. We apply three different approaches to generate the matrix S: (1) Uniformly distributed (UD) sensor placement, which requires no prior knowledge of the flow and places all the sensors uniformly across the observation window; (2) Modeluncertainty-based (MU) sensor placement, which places sensors in regions where the largest model uncertainty occurs without filtering; and (3) pivoted QR sensor placement [18], which optimally chooses the sensor locations based on the dominant POD modes obtained from the slow measurements. More details about these approaches and the sensor locations are described in Appendix A. We place 16 sensors in the simulations and compare the reconstructions as a representative case to evaluate these three placement strategies. In addition, we compare the reconstruction errors of QR placement with 9 sensors vs. 16 sensors to study the influence the number of the sensors has on the reconstruction. The fast measurements are assumed to be temporally resolved with the same sampling time as the DNS snapshots, i.e., ∆t + f = t + = 0.065. We collect the data from JHTDB with a total time duration T + = 37.96, yielding a total of 584 snapshots. We applied singular value decomposition to the first 200 snapshots to tailor the POD basis for QR pivoting approach in sensor placement as described in Appendix A, and the remaining 384 snapshots are subsampled to simulate the measurements for reconstruction, as well as serving as ground truth for validation of the fusion algorithm. This results in four slow-time intervals for investigating the reconstruction.

Reconstruction error evaluation
The Root Mean Square Error (RMSE) of the reconstruction can be calculated as whereû 1 andû 2 are the reconstructed velocity fluctuations while u 1 and u 2 are the velocity fluctuation from the DNS ground truth. Recall that x 2 = 0 indicates the location at the lower wall of the channel, and x 2 = h represents the centerline of the channel. The reconstruction accuracy is evaluated for forward estimates, backward estimates, and the forward-backward estimates using different weighting schemes as described in Section 2.2.
To analyze the impact of uncertainty from the measurement and process noise, we conduct Monte Carlo simulations and evaluate the mean and variance of the RMSE over multiple realizations. Following the same metrics, we can also evaluate the RMSE with respect to the streamwise and wall-normal velocity components separately, as follows: The equations above aim to evaluate the reconstruction errors for each component of the 2D velocity profiles. Specific details will be presented with the results in Section 4.

Results and Discussion
In this section, we evaluate the reconstruction performance of the proposed fusion algorithms and investigate the influence of implementation decisions (e.g., sensor placement, weighting functions). Section 4.1 compares the time evolution of the RMSE of the forward, backward, and fused estimates. Section 4.2 compares these reconstructed flow fields to model prediction without sensor fusion and to the JHTDB ground truth. Section 4.3 evaluates the influence of the number and placement of sensors on reconstruction accuracy. In Section 4.4, the reconstructed turbulence intensity and Reynolds stresses are compared to DNS to further evaluate the fusion performance. Unless otherwise state, the total number of sensors will be 16 for all results reported. Corresponding sensor locations are reported in the Appendix. Table 1 summarizes the different turbulent flow reconstruction techniques as well as the sensor placement approaches used in this study. Figures with the corresponding results are also listed in Table 1.

Flow reconstruction RMSE comparison
We begin by assessing the reconstruction performance of the various fusion algorithms using Monte Carlo simulations. The error statistics of the temporally weighted smoother (KS t ± ) and spatiotemporally weighted smoother (KS st ± ) are compared with those of the forward filter (KF + ) and the backward filter (KF − ). The mean and 1-σ (1 standard deviation) bounds from 18 independent realizations over four slow sampling intervals (∆t + s = 6.24) are reported in Fig. 5. These results correspond to the case of 16 fast sensors arranged according to the QR placement approach. An initial RMSE of ≈ 0.05 is observed at the beginning and end of each interval, and can be at- tributed to the noise floor determined by the sensor noise that contaminates the slow snapshot measurements. The maximum RMSE of the fused forward-backward reconstruction is found to be ≈ 0.2 over each interval, regardless of the weighting scheme employed.
Comparing the forward estimate with the backward estimate, we observe that the backward prediction tends to exhibit a larger RMSE than the forward prediction over each interval. The same observation was made regarding backward predictions from the RDT model without fast sensors or fast filtering, as reported in our previous study [17]. As noted earlier, this can be explained by the fact that the backward RDT model contains a negative viscous diffusion term which leads to a steep increase in prediction error over time. As expected, the error in the forward filter increases toward the end of an interval, while the error in the backward filter increases toward the beginning of an interval. The smoothing algorithms combine these forward and backward estimates in a manner that weighs the more reliable estimate-at a particular time instant or point in space-more heavily, thus yielding an improvement in the reconstruction performance over the entirety of each interval.
From Fig. 5, it can also be observed that the reconstruction accuracy varies between the different intervals. These differences are especially prominent in the estimates from the forward filter (KF + ). This indicates that filter initialization plays a role in reconstruction performance, as is to be expected. We plot the mean and 1-σ bounds of interval-averaged RMSE in Fig. 6, based on the 18 realizations over the four intervals (i.e., 72 total realizations). Accounting for the reconstruction performance in this way reveals a wider band in the estimation error variance. However, the magnitude of the RMSE never exceeds ≈ 0.2 for reconstructions based on fusion of the forward and backward estimates. This suggests that reconstruction performance is not restricted to the specific dataset or flow profile used in this study. Another interesting observation is that smoothing with the temporal weighting scheme (KS t ± ) performs slightly better than smoothing with the spatiotemporal weighting scheme (KS st ± ), at least for the time-horizon and data considered here. This can be seen in Fig. 6, where the temporal weighting scheme generates predictions with a narrower variance than the spatiotemporal weighting scheme. In addition, the spatiotemporal weighting is also observed to exhibit a large error at t + ≈ 0.5. These performance differences can be attributed to the fact that the spatiotemporal weighting scheme is based purely on the flow physics, not taking into consideration any specific properties concerning the specific configuration of point sensors. Accounting for the placement of fast point sensors could improve performance here; however, obtaining and implementing such a weighting scheme would be substantially more complicated compared to a simple weighting scheme based on the ROI/DOD figures shown in Fig. 4. This is precisely the reason that the spatiotemporal weighting scheme works best when applied to the uniformly distributed (UD) sensor placement, as we will see later in Section 4.3. We can consider the temporal and spatiotemporal weighting schemes to be more "universal" in the context of multi-sensor fusion, since they are agnostic to the spatial positioning of point sensors in the domain.
The plots of x1 (t) and x2 (t) in the time domain are reported in Fig. 7. The fact that these reconstruction errors are normalized with respect to the spatial integration of independent velocity components (i.e., u 1 and u 2 ) leads to different RMSE magnitudes between the two. The maximum streamwise RMSE for both smoothers is approximately x1 ≈ 0.15, while the maximum wall-normal RMSE for both smoothers is approximately x2 ≈ 0.3. We note that the reconstruction of the wall-normal velocity fluctuations u 2 is more sensitive to measurement noise than the reconstruction of the streamwise velocity fluctuations u 1 . This sensitivity to measurement uncertainty arises on account of the artificial viscosity that is introduced from discretizing the governing equations via finite differences. This artificial viscosity serves to dampen higher-frequency and smaller-scale components more than lower-frequencies and larger-scale components. For wall-normal velocity u 2 , higher energy appears at the higher frequencies than for the streamwise velocity u 1 ; thus, the attenuation is larger for wall-normal velocity u 2 . This explains the fact that the wall-normal RMSE x2 (t) is found to be larger than the streamwise RMSE x1 (t). Up to this point, we have posited that multi-sensor fusion can improve the accuracy of the reconstructed flow compared with reconstructions from model predictions based on non-time-resolved field measurements alone. To see that this is true, we compare the resulting RMSE for the multi-sensor fusion approach with those from model predictions, shown in in Fig. 8. The multi-sensor fusion approachregardless of weighting scheme-yields a lower mean RMSE compared to the modelprediction-based reconstructions. As an example, consider the third interval (i.e., 12.48 ≤ t + ≤ 18.72), the RMSE from temporal weighting of the forward-backward model predictions (RDT t ± ) reaches approximately ≈ 0.27, whereas multi-sensor fusion with the same temporal weighting (KS t ± ) yields a maximum RMSE of ≈ 0.17-an improvement of about 40%. Even for the other intervals, the percentage improvement is approximately 20%. However, comparing results between the spatiotemporal weighting schemes is less encouraging. Here, the performance gains achieved by multi-sensor fusion with spatiotemporal weighting (KS st ± ) are less prominent, but do outperform the spatiotemporal model prediction (RDT st ± ) in terms of the mean RMSE in each interval. For model predictions, spatiotemporal weighting yields a lower maximum error than temporal weighting. This is in contrast to the multisensor fusion results, for which the temporal and spatiotemporal weighting scheme does not have noticeable effect on reconstruction accuracy. As discussed above, the spatiotemporal weighting scheme was not designed with the spatial locations of the fast measurements in mind. It may be possible to achieve superior reconstruction ac- curacy with an alternative weighting scheme that takes the sensor arrangement into account.
Considering the statistics between all realizations and all intervals in a single Monte Carlo plot better highlights the aggregate performance between these approaches (see Fig. 8(b)). An interesting observation here is that the RMSE of the model-prediction-based reconstruction tends to possess a smaller variance compared to the multi-sensor fusion result. The model-prediction-based reconstructions only make use of noisy snapshot data at the beginning and end of each interval. In contrast, the multi-sensor fusion approach makes use of these same noisy snapshots, but also noisy measurements from the fast sensors. Although the fast measurements provide time resolved information for the fast filter, the sensor noise from these measurements introduces uncertainty into the estimate at each fast time-step. This uncertainty results in a larger estimation variance, which compounds over time.

Flow field reconstruction comparison
In the previous section, we reported on the error statistics of the flow reconstruction based on a Monte Carlo analysis over multiple realizations. Now, we consider the flow field reconstruction for a single representative realization more closely. We begin by comparing performance between the forward filter (KF + ) and the two smoothers, one with the temporal weighting scheme (KS t ± ) and the other with the spatiotemporal weighting scheme (KS st ± ). Figs. 9(a)-(b) show the evolution of the ground truth and reconstructed velocity fluctuations along a wall-normal slice of the channel at a streamwise station of x 1 = 0.5h. These results correspond to four consecutive slow-sampling intervals and the case of 16 fast sensors, arranged according to the QR placement method described in Appendix A.
From visual inspection, it appears that both smoothers perform comparably, and capture many of the large and small scale features appearing in the ground truth. The forward filter also captures dominant features in the flow, but finer details are not captured as well. To better quantify and assess this performance, we also report the reconstruction errors in Figs. 9(c)-(d). These are obtained by subtracting the ground truth from each reconstruction. As discussed when assessing performance in terms of RMSE, multi-sensor smoothing with either temporal or spatiotemporal weighting schemes perform comparably in a Monte Carlo study of RMSE. Here too, we find that these two weighting schemes yield reconstructions with comparable error performance in time and space. Both smoothers outperform the forward filter in terms of overall error, but the dominant errors appear to be concentrated at the same points in space-time, at least for this slice of the flow.
Note that the reconstruction errors for the fast filter tend to be smaller immediately after a slow snapshot is made available, then grow in time up until the next snapshot. From Fig. 9, it can be observed that flow field reconstruction from the forward filter is less smooth towards the end of each slow sampling interval, which leads to slightly larger reconstruction errors at the end of the interval. This observation is consistent with the relatively large RMSE at the end of each interval, as reported in Fig. 5. The flow field reconstruction from the smoothers do not exhibit the same degree of nonsmoothness toward the end of the interval: reconstruction error grows at first, but then decreases with time until the end of the interval. These differences with the forward filter and the two smoothers are to be expected, since the smoother was designed to improve upon the forward filter results by taking a backward pass through the data using a backward filter.
The results above compare flow reconstructions obtained by fusing measurements from slow and fast sensors. The associated reconstruction performance can be benchmarked against the case of reconstructing the flow from model predictions based on only the slow measurements. Figs. 10(a)-(b) report the flow reconstruction over the same wall-normal slice of the channel as before, but now showing results from spatiotemporally weighted forward and backward model predictions (RDT st ± )-i.e., the model-prediction approach with the best reconstruction performance, based on the RMSE results reported in Fig. 8. Results for the spatiotemporally weighted multisensor smoother (KS st ± ) and the ground truth are also reported here. The reconstruction based on the multi-sensor fusion is notably closer to the ground truth than the  model-prediction-based reconstruction. These differences are more clearly seen in a visualization of the associated reconstruction errors shown in Figs. 10(c)-(d). From this analysis, it is evident that multi-sensor fusion yields a lower reconstruction error than the model prediction approach, both in the streamwise and the wall-normal components of the velocity fluctuations. Interestingly, even though the model-predictionbased reconstructions are based on both a forward and backward pass on the data, this approach tends to exhibit larger errors closer to the starts and ends of the slow sampling intervals than the multi-sensor approach. Indeed, the fast sensors provide valuable information that benefits the reconstruction across the entire interval, including at times that the slow snapshots are available.
Finally, we investigate the reconstructions from RDT st ± and KS st ± as a function of both spatial directions, focusing on a single time instant t + = 3.12 that is halfway between two slow snapshots (see Fig. 11). Note that this is the time instant within an interval corresponding to relatively large RMSE in all of the reconstruction approaches (see Fig. 5). The most notable difference arises in the reconstruction of wall normal velocity fluctuations. The multi-sensor fusion approach undoubtedly captures the dominant features and many of the finer features that can be seen in the ground truth, whereas the model-prediction-based reconstruction smears out many of these features. The model-prediction approach performs better for reconstructing the streamwise velocity fluctuations than for wall-normal velocity fluctuations. However, the multi-sensor fusion approach still captures finer features that the model prediction approach misses. These observations are consistent with the fact that wall-normal re-construction error x2 was found to be greater than the streamwise reconstruction error x1 . Further, these results suggest that the dominant contributor to the RMSE improvements realized by the multi-sensor fusion approach can be attributed to improvements in reconstructing the wall-normal velocity fluctuations. The fast sensors provide sufficient information that can be exploited by the fusion framework to improve model predictions based on RDT st ± from a maximum mean RMSE of ≈ 0.23 down to ≈ 0.17 for KS st ± (see Fig. 8).

Effect of number and arrangement of fast sensors
All of the multi-sensor fusion results reported to this point have corresponded to the case of m = 16 fast sensors distributed in the channel based on the pivoted QR sensor arrangement described in Appendix A. Here, we briefly investigate the influence of number of fast sensors m and their arrangement on the quality of the flow reconstruction. We consider sensor arrangements determined using three different approaches, all described in Appendix A: (1) uniformly distributed (UD) sensor placement, (2) model-uncertainty-based (MU) sensor placement, and (3) pivoted QR sensor placement. We only consider perfect squares m = 9 and m = 16 here in order to satisfy symmetry conditions for the UD placement. Both values of m are roughly three orders of magnitude less than the dimension of the flow state being reconstructed. In the course of our study, we encountered difficulties in filter tuning when m < 9 fast sensors were utilized, and so do not report on those results here. Although reconstruction performance can be further improved through the use of a larger number of sensors (i.e., m > 16), our ultimate interest is to leverage multisensor fusion for flow reconstruction in actual experiments where a larger number of fast sensors may be impractical and so do not report on such cases here either. In Fig. 12, we report the RMSE for reconstruction using UD, MU, and QR sensor arrangements corresponding to m = 9 and m = 16 fast sensors. Results are reported based on a Monte Carlo study with 72 realizations. The mean RMSE from this Monte Carlo study is plotted for each case using a dotted line, and the shaded regions correspond to the 1-σ (i.e., one standard deviation) bounds in the reconstruction resulting from the 72 realizations. From these results, it is evident that a larger number of sensors results in a lower RMSE, regardless of the sensor arrangement-i.e., m = 16 sensors yields a lower RSME on average than m = 9 sensors for a given sensor placement strategy. Similarly, a larger number of sensors tends to yield a lower variance as well. This is especially evident in the QR placement results. Our general finding that more sensors improves performance, both in terms of reducing RMSE and reducing variance, is to be expected. However, it is interesting to note that the reconstruction variance is found to be least sensitive to the number of sensors when the UD sensor placement strategy is employed. This indicates the UD sensor arrangement is less sensitive to the flow conditions and noise than the other sensor arrangements. We also note that for both the UD and MU sensor arrangements, a spatiotemporal weighting scheme within the smoother improves the RMSE and variance of the reconstruction relative to a temporal weighting scheme. In contrast, in the case of the QR sensor arrangements, both the RMSE and variance of the reconstruction tend to be lower when a temporal weighting scheme is used, as opposed to a spatiotemporal weighting scheme, in the smoothing algorithm; however, this performance difference is less pronounced in the case of m = 16 than it is for m = 9. The results here suggest that future work on determining an optimal sensor arrangement should consider the balance between RMSE and variance in the reconstruction error, and the sensitivities of each to specific sources of uncertainty.
The results of the Monte Carlo study described here indicate that sensor fusion based on the QR sensor arrangement will yield the best reconstruction performance with a temporal weighting scheme in the smoother; whereas, for UD and MU sensor arrangements, a spatiotemporal weighting scheme will yield the better performance. In Fig. 13(a)-(b), we report representative results for the reconstructed velocityfluctuation-field based on each of these three sensor arrangements with m = 16 and the corresponding weighting scheme yielding the best reconstruction performance. It is evident that the reconstructed flow field from all three approaches captures the dominant flow features observed in the JHTDB ground truth. The QR sensor arrangement appears to yield a better match to the ground truth than either the UD or MU arrangements, especially for the streamwise component of velocity. Reconstructions from both the UD and MU arrangements appear strikingly similar to each other. To better quantify the performance, we report the associated reconstruction error relative to the JHTDB ground truth in Fig. 13(c)-(d). These results confirm that the QR sensor arrangement outperforms the other two sensor arrangements by yielding a lower magnitude of reconstruction error. To explain this difference, recall that the QR placement procedure (see Appendix A) made use of a tailored basis of POD modes-extracted from 200 snapshots of the flow field-to determine the optimal arrangement of fast point sensors. Although the snapshots used to determine the POD basis are distinct from the snapshots used to evaluate the sensor fusion algorithms, it still stands that the tailored basis of POD modes allows the QR pivoting procedure to determine the sensor locations that will best capture fluctuations in the turbulent kinetic energy. Indeed, in the next section we will investigate this idea further, and show that the QR sensor arrangement also outperform the other placement strategies in terms of capturing the second order statistics of the turbulent fluctuations. In contrast, the UD and MU placement strategies do not consider the statistical characteristics of the turbulent fluctuations. The UD strategy takes no physics into account at all. The MU method only considers placement of sensors into high-error regions of the flow based on the model-based reconstruction. The variance of these errors was not considered for the MU placement, but may be a useful heuristic to consider for fast point sensor placement in future investigations.

Second order statistics comparison
We now assess reconstruction performance of multi-sensor fusion in terms of capabilities for reproducing the second order statistics of the JHTDB ground truth flow. In the previous section, we found that the QR arrangement of fast sensors modestly outperformed the UD and MU arrangements in terms of RMSE from a Monte Carlo study. We expect that multi-sensor fusion based on the QR sensor arrangement will out-xxiv Mengying Wang et al.  perform reconstructions of the second order flow statistics based on UD and MU arrangement. This is because the approach for determining the QR arrangement makes use of energy-optimal modes from the proper orthogonal decomposition (POD). Indeed, this turns out to be the case, as can be seen in Fig. 14.
Fig. 14: Comparison of reconstructed second-order statistics u 2 1 (left), u 2 2 (middle) and u 1 u 2 (right) spatially averaged along the wall-normal direction and over 72 realizations. Note that the statistics do not resemble those for canonical channel flow due to the short time period being considered, i.e., they are not fully converged. Fig. 14 shows the turbulence intensities u 2 1 and u 2 2 along with the Reynolds stress u 1 u 2 as a function of wall-normal station x + 2 for the JHTDB ground truth and reconstructions based on model predictions (RDT t ± , RDT st ± ) and multi-sensor fusion using on UD, MU, and QR sensor arrangements (UD KS st ± , MU KS st ± , QR KS t ± ). Far from the wall (approximately x + 2 > 100), all methods capture the second order statistics comparably well. Closer to the walls, the deviation between the reconstructed statistics and the ground truth is found to be larger. This indicates that strong turbulence near the walls tends to degrade the reconstruction performance. However, there is a notable difference in overall performance between multi-sensor fusion using the QR sensor arrangement and the other methods. The temporally weighted smoother with QR sensor placement (dashed red line) demonstrates a markedly superior ability to capture the turbulence intensities u 2 1 and u 2 2 of the ground truth relative to the other reconstruction approaches for all x + 2 -exhibiting about half the deviation from the ground truth as the other approaches. This is also true for reconstruction of the Reynolds stress u 1 u 2 ; however, the ability to recover the Reynolds stress near the wall using the QR arrangement tends to be less pronounced. In contrast, multisensor fusion with UD fast sensor arrangement yields comparable reconstruction of the second-order statistics as the model-prediction-based methods.
Considering the second-order statistical analysis of these three sensor arrangements with the corresponding weighting scheme yielding the best reconstruction performance, QR sensor placement with the temporal weighting scheme outperforms the other two placements with spatiotemporal weighting scheme, especially in reproducing the turbulence intensities. These results highlight the utility and related energy-optimal sensor placement strategies for turbulent flow reconstruction, even within the context of multi-sensor fusion.

Conclusion
In this paper, we have proposed and evaluated methods for model-based multi-sensor fusion for turbulent flow reconstruction. We showed that non-time-resolved planar field measurements, time-resolved point measurements, and predictions from a simple physics-based model can be combined to reconstruct the velocity fluctuations in a turbulent channel flow from the Johns Hopkins Turbulence Database (JHTDB). This was achieved by using a "fast filter" to fuse fast point sensor measurements with RDT-based model predictions, then to use a "slow filter" to fuse these fast filter estimates with the non-time-resolved field measurements on a slower time-scale. We further showed that the multi-rate filtering scheme can be used both forward and backward in time, and that physically motivated weighting schemes can be used to yield a smoothing algorithm with reduced error and variance in the flow reconstruction. In addition, we found that these smoothing algorithms were also capable of reproducing the second order statistics of the underlying turbulent flow.
Monte Carlo simulations were used to investigate the statistical performance of the multi-sensor fusion methods. The influence of the number and arrangement of fast point sensors was also investigated, with three sensor placement approaches studied: (1) uniformly distributed (UD) sensor placement, (2) model-uncertainty-based (MU) sensor placement, and (3) pivoted QR sensor placement. For all methods, reproduction quality improved with a larger number of sensors. However, among the three methods, the QR sensor placement strategy was found to have the best performance overall, yielding a lower RMSE and a better reproduction of the turbulence intensities. The UD scheme is easier to implement a priori and could also be applied in reconstruction with multi-rate PIV systems, e.g., high-rate low-resolution PIV and low-rate high-resolution PIV.
The results of our study draw attention to several avenues worthy of future research and investigation, especially with regards to enabling multi-sensor fusion for turbulent flow reconstruction in practice. Firstly, the need for computationally efficient multi-sensor fusion algorithms cannot be understated. The propagation of the covariance matrix within the slow and fast filters consumes the bulk of the computational time and effort. The covariance matrix is of the order N 2 , where N is the state dimension. Although the covariance matrix is symmetric-a structural property that can be exploited for computational speed-up-it stands that this will be the limiting computational bottleneck in larger-scale problems. In addition, sensor noise models used in this study were relatively simplified. Future investigations must consider the role of sensor noise that are representative of the specific instrumentation being employed for flow reconstruction. For example, the measurement uncertainty associated with the velocity field measured using a PIV system will necessarily be a function of the velocity magnitude. Thus, simple additive Gaussian noise models for sensor noise do not account for these more realistic uncertainty profiles-even though they are standard choices for algorithm development in the literature on flow reconstruc-tion. Finally, the role of the number and arrangement of fast sensors is central to achieving successful flow reconstruction. As we have seen in this study, some arrangements are better suited for reducing RMSE and improving the reconstruction of second-order statistics (e.g., pivoted QR placement); however, we have also seen that some arrangements yield lower variance in the flow reconstruction and thus exhibit less sensitivity/greater robustness to operating conditions (e.g., uniformly distributed sensors). Further the weighting schemes used for smoothing between the forward and backward filter estimates did not directly account for the placement of sensors. Future work must consider the placement of sensors and the weighting functions used for smoothing simultaneously. We expect that by doing so, multi-sensor fusion will yield great improvements in flow reconstruction performance, in terms of minimizing reconstruction errors and variances, but also in terms of capturing a wider range of turbulent physics.

Appendix A Sensor placement
In this study, three different approaches are investigated for determining the arrangement of fast point sensors-encoded in the term S-for multi-sensor fusion: (1) uniformly distributed (UD) sensor placement, (2) model-uncertainty-based (MU) sensor placement, and (3) pivoted QR sensor placement. Each of these is described here.
The UD placement approach is the simplest of the three approaches considered. UD placement yields a uniform distribution of m fast point sensors throughout the domain. Sensors are evenly spaced from each other and from the domain boundaries along both the x 1 and x 2 directions. Although this approach does not take into account any flow physics, it provides a reasonable benchmark for comparison. The same UD approach was also considered in [30]. In order to maintain symmetry in the resulting placement, the number of sensors m is selected to be a perfect square. Here, we consider only m = 9 and m = 16, which are values determined for by the QR placement approach-to be described momentarily-and used here for direct comparison. The UD sensor placements for these two cases are plotted as black stars in Fig. 15, overlaid on a snapshot of the streamwise velocity fluctuations u 1 at an arbitrary instant in time.
The UD method implicitly assumes that velocity fluctuations at all spatial points are equally important for flow reconstruction, which is not true in general. A simple means of addressing this shortcoming is to consider placing fast point sensors at locations with high reconstruction error, when a model-based reconstruction is used directly-i.e., a model-uncertainty-based (MU) placement strategy. In Fig. 16, we report wall-normal variation in error over time (x 2 , t) between the RDT-based reconstruction and JHTDB ground truth averaged over 50 realizations of the flow. The wall-normal variation in error over time is defined as  The final approach considered here is the sparse sensor placement algorithm based on a pivoted QR decomposition, as presented in [18]. The QR approach uses  snapshots of the flow field to obtain a tailored (reduced) basis of POD modes for capturing dominant signals in the flow. Then, a column-pivoted QR factorization is used to construct S, with the leading m column-pivots indicating the m sensor locations that will best approximate the training data in a least-squares sense. The specific formulation and additional details of the approach can be found in [18]. We emphasize that this QR approach requires training data for computing a tailored basis of POD modes, but that temporally resolved velocity field data will not be available in practice: this was the motivation for multi-sensor fusion in the first place. As such, we collect 200 slow-in-time snapshots of the flow field to extract a POD basis. The marginal contribution of each new POD mode to the energy in the training data is reported in Fig. 18. This cumulative energy analysis indicates that a tailored basis of r = 16 POD modes captures roughly 83% of the energy. Thus, a QR-based sensor placement on this tailored basis will require m = 16 sensors. We also consider a tailored basis of r = 9 POD modes (captured using m = 9 sensors) for comparison with the UD approach, which requires the number of sensors to be a perfect square to maintain symmetry. We note that the specific sensor arrangements resulting from the QR approach were found to be sensitive to the training data. This sensitivity was observed to be more significant with a larger number of sensors. However, even with sensitivities and variations in the sensor arrangements, a larger number of sensors was found to improve flow reconstruction so long as appropriate filter tuning was performed within the multi-sensor fusion framework. Due to these additional complexities, we only consider the cases of m = 9 and m = 16 sensors in this study.
The sparse sensor arrangements resulting from the subsequent column-pivoted QR procedure for m = 9 and m = 16 sensors are shown in Fig. 19. Most of the sensors in these arrangements are placed in the vicinity of the lower wall. This confirms that the streamwise velocity fluctuations are a dominant feature of this flow Normalized cumulative energy 0.83 Fig. 18: Cumulative contribution of the leading r POD modes to the total energy of the training data. The leading 16 POD modes capture 83% of the energy in the training data, and the leading 9 POD modes capture 68%. that needs to be captured for flow reconstruction. Interestingly, near wall information was also found to be important using the MU placement approach. However, the MU approach resulted in a uniform streamwise distribution of sensors at a wall-normal station of x + 2 ≈ 100; in contrast, the QR placement tends to yield sensor locations with x + 2 < 100 and an irregular streamwise distribution of sensors, which would not have been obtained using simple heuristics.  Friction velocity v f ∼ N (0, R f ) Point measurement noise with Gaussian distribution of zero mean and covariance matrix R f v s ∼ N (0, R s ) Field measurement noise with Gaussian distribution of zero mean and covariance matrix R s w ∼ N (0, Q) Process noise with Gaussian distribution of zero mean and covariance matrix Q (x 1 , x 2 , x 3 ) Three dimensional coordinate system y f "Fast" point measurements y s "Slow" field measurements