2.1. Dimensionality reduction

We use the snapshot method of Sirovich (Reference Sirovich1987) to obtain the dominant POD modes for the velocity field. The snapshot matrix $ \boldsymbol{Y}\left(\boldsymbol{x},{t}_1:{t}_K\right) $ for the velocity fluctuations $ {u}^{\prime } $ and $ {v}^{\prime } $ is assembled:

(2.1) $$ \boldsymbol{Y}\left(\boldsymbol{x},{t}_1:{t}_K\right)=\left[\begin{array}{cccc}{u}^{\prime}\left({\boldsymbol{x}}_1,{t}_1\right)& {u}^{\prime}\left({\boldsymbol{x}}_1,{t}_2\right)& \dots & {u}^{\prime}\left({\boldsymbol{x}}_1,{t}_K\right)\\ {}\vdots & \vdots & \dots & \vdots \\ {}{u}^{\prime}\left({\boldsymbol{x}}_N,{t}_1\right)& {u}^{\prime}\left({\boldsymbol{x}}_N,{t}_2\right)& \dots & {u}^{\prime}\left({\boldsymbol{x}}_N,{t}_K\right)\\ {}{v}^{\prime}\left({\boldsymbol{x}}_1,{t}_1\right)& {v}^{\prime}\left({\boldsymbol{x}}_1,{t}_2\right)& \dots & {v}^{\prime}\left({\boldsymbol{x}}_1,{t}_K\right)\\ {}\vdots & \vdots & \dots & \vdots \\ {}{v}^{\prime}\left({\boldsymbol{x}}_N,{t}_1\right)& {v}^{\prime}\left({\boldsymbol{x}}_N,{t}_2\right)& \dots & {v}^{\prime}\left({\boldsymbol{x}}_{N,{t}_K}\right)\end{array}\right], $$

where $ \boldsymbol{Y}\in {\mathrm{\mathbb{R}}}^{2N\times K} $ , $ {\boldsymbol{x}}_i=\left[{x}_i,{y}_i\right] $ $ \left(i=1,2,\dots N\right) $ is the location vector for the $ i $ -th spatial location (cell centroid), $ N $ is the number of cells, and $ K $ the number of snapshots (we assume $ K<2N $ ). Singular value decomposition (SVD) is performed on the weighted matrix $ {\mathcal{V}}^{1/2}\boldsymbol{Y} $ ,

(2.2) $$ {\mathcal{V}}^{1/2}\boldsymbol{Y}={\boldsymbol{U}}_Y{\varSigma}_Y{\boldsymbol{V}}_Y^{\top }, $$

where $ \mathcal{V}=\operatorname{diag}\left({V}_1,{V}_2\dots {V}_N,{V}_1,{V}_2,\dots, {V}_N\right) $ is a diagonal matrix with the cell volumes $ {V}_i $ in the main diagonal, $ {\boldsymbol{U}}_Y\in {\mathrm{\mathbb{R}}}^{2N\times K} $ contains the left singular vectors, $ {\varSigma}_Y\in {\mathrm{\mathbb{R}}}^{K\times K} $ is a diagonal matrix that stores the singular values, and $ {\boldsymbol{V}}_Y\in {\mathrm{\mathbb{R}}}^{K\times K} $ contains the right singular vectors. The scaled POD eigenmodes $ {U}_{Y,k}\left(\boldsymbol{x}\right) $ are extracted from the columns of $ {\boldsymbol{U}}_Y $ :

(2.3) $$ {U}_Y\left(\boldsymbol{x}\right)={\mathcal{V}}^{-1/2}{\boldsymbol{U}}_Y, $$

and the time coefficients are extracted from

(2.4) $$ \boldsymbol{a}(t)={\boldsymbol{\Sigma}}_Y{\boldsymbol{V}}_Y^{\top }. $$

The fluctuating velocity field can be written as

(2.5) $$ {u}_i^{\prime}\left(x,y,t\right)=\sum \limits_{k=1}^K{a}_k(t){U}_{Y,k}^{(i)}\left(x,y\right)\approx \sum \limits_{k=1}^{m_u}{a}_k(t){U}_{Y,k}^{(i)}\left(x,y\right),\left(i=1,2\right), $$

where $ {m}_u\left(\ll K\right) $ is the number of retained POD modes, $ {a}_k(t) $ is the time coefficient of the $ k $ -th POD mode, and $ {U}_{Y,k}^{(i)} $ is the $ k $ -th POD eigenmode of the $ i $ -th velocity component. This expression can be written in the matrix form as

(2.6) $$ \underset{\equiv {\boldsymbol{u}}^{'}\left(\boldsymbol{x},t\right)}{\underbrace{\left[\begin{array}{l}{u}_1^{\prime}\\ {}{u}_2^{\prime}\end{array}\right]}}=\underset{\equiv {\boldsymbol{U}}_Y\left(\boldsymbol{x}\right)}{\underbrace{\left[\begin{array}{cccc}{U}_{Y,1}^{(1)}& {U}_{Y,2}^{(1)}& \dots & {U}_{Y,{m}_u}^{(1)}\\ {}{U}_{Y,1}^{(2)}& {U}_{Y,2}^{(2)}& \dots & {U}_{Y,{m}_u}^{(2)}\end{array}\right]}}\underset{\equiv \boldsymbol{a}(t)}{\underbrace{\left[\begin{array}{c}{a}_1\\ {}{a}_2\\ {}\vdots \\ {}{a}_{m_u}\end{array}\right]}} $$

or more compactly

(2.7) $$ {\boldsymbol{u}}^{\prime}\left(x,y,t\right)={\boldsymbol{U}}_Y\left(x,y\right)\boldsymbol{a}(t). $$

2.2. Computation of the time-delayed Hankel coefficients

Using the POD coefficients from the training dataset, we assemble the time-delayed Hankel matrix $ \boldsymbol{H} $ :

(2.8) $$ \boldsymbol{H}=\left[\begin{array}{cccc}\boldsymbol{a}\left({t}_1\right)& \boldsymbol{a}\left({t}_2\right)& \dots & \boldsymbol{a}\left({t}_p\right)\\ {}\boldsymbol{a}\left({t}_2\right)& \boldsymbol{a}\left({t}_3\right)& \dots & \boldsymbol{a}\left({t}_{p+1}\right)\\ {}\vdots & \vdots & \ddots & \vdots \\ {}\boldsymbol{a}\left({t}_q\right)& \boldsymbol{a}\left({t}_{q+1}\right)& \dots & \boldsymbol{a}\left({t}_{K_{\mathrm{train}}}\right)\end{array}\right], $$

where we use $ q $ vectors of $ {m}_u $ POD coefficients in each column. The number of columns is $ p={K}_{\mathrm{train}}-q+1 $ , where $ {K}_{\mathrm{train}} $ is the number of snapshots for the training dataset; thus, $ \boldsymbol{H}\in {\mathrm{\mathbb{R}}}^{\left({m}_u\times q\right)\times p} $ . Performing SVD on $ \boldsymbol{H} $ we obtain,

(2.9) $$ \boldsymbol{H}\approx {\boldsymbol{U}}_H{\boldsymbol{\Sigma}}_H{\boldsymbol{V}}_H^{\top }, $$

where $ {\boldsymbol{U}}_H\in {\mathrm{\mathbb{R}}}^{\left({m}_u\times q\right)\times r} $ , $ {\boldsymbol{\Sigma}}_H\in {\mathrm{\mathbb{R}}}^{r\times r} $ , $ {\boldsymbol{V}}_H\in {\mathrm{\mathbb{R}}}^{p\times r} $ and $ r $ is the number of retained singular values. The matrix $ {\boldsymbol{U}}_H $ of the left singular vectors can be explicitly written as

(2.10) $$ {\boldsymbol{U}}_H=\left[\begin{array}{cccc}{\boldsymbol{U}}_{H,1}^{\left(u,v\right)}\left({t}_1\right)& {\boldsymbol{U}}_{H,2}^{\left(u,v\right)}\left({t}_1\right)& \dots & {\boldsymbol{U}}_{H,r}^{\left(u,v\right)}\left({t}_1\right)\\ {}{\boldsymbol{U}}_{H,1}^{\left(u,v\right)}\left({t}_2\right)& {\boldsymbol{U}}_{H,2}^{\left(u,v\right)}\left({t}_2\right)& \dots & {\boldsymbol{U}}_{H,r}^{\left(u,v\right)}\left({t}_2\right)\\ {}\vdots & \vdots & \ddots & \vdots \\ {}{\boldsymbol{U}}_{H,1}^{\left(u,v\right)}\left({t}_q\right)& {\boldsymbol{U}}_{H,2}^{\left(u,v\right)}\left({t}_q\right)& \dots & {\boldsymbol{U}}_{H,r}^{\left(u,v\right)}\left({t}_q\right)\end{array}\right], $$

where each column $ {\boldsymbol{U}}_{H,i}^{\left(u,v\right)}\left({t}_1:{t}_q\right)\in {\mathrm{\mathbb{R}}}^{m_u} $ is the $ i $ -th delayed singular mode. The matrix of the right singular vectors $ {\boldsymbol{V}}_H $ can be explicitly written as

(2.11) $$ {\boldsymbol{V}}_H=\left[\begin{array}{cccc}{v}_{H,1}\left({t}_1\right)& {v}_{H,1}\left({t}_2\right)& \dots & {v}_{H,1}\left({t}_p\right)\\ {}{v}_{H,2}\left({t}_1\right)& {v}_{H,2}\left({t}_2\right)& \dots & {v}_{H,2}\left({t}_p\right)\\ {}\vdots & \vdots & \ddots & \vdots \\ {}{v}_{H,r}\left({t}_1\right)& {v}_{H,r}\left({t}_2\right)& \dots & {v}_{H,r}\left({t}_p\right)\end{array}\right]. $$

2.3. Construction of a dynamical system for POD coefficients

The next step is the construction of a linear dynamical system that describes the evolution of the POD coefficients. We first describe two algorithms, linear 1 and linear 2, that extract linear models for the right singular vectors in matrix $ {\boldsymbol{V}}_H $ . A schematic of the two algorithms is presented in Figure 1. We then briefly describe an algorithm based on the n4sid system identification method.

Figure 1.

Schematic overview of the algorithms linear 1 and linear 2.

Figure 1. Long description

The diagram is organized into several stages of data processing.

1. Top Section: Flow field data is shown as a stack of three heatmaps labeled k equals 1, k equals 2, and k equals K, representing Y(x, t sub k). An arrow points to an Economy S V D section where V super 1/2 Y is approximately equal to U sub Y, Sigma sub Y, and V sub Y super T. Dimensions are labeled as 2 N, K, and m sub u.

2. Middle Left Section: True output a[k] is shown as four oscillating waveforms for m equals 1 through 4. Below this, four heatmaps show u prime and v prime components for U sub Y,1 and U sub Y,2. At the bottom left, sensor measurements for u prime and v prime are plotted as red and blue waves, labeled as sensor input s[k].

3. Middle Right Section: The Economy S V D output feeds into Delayed coordinates, forming a matrix H. This matrix undergoes another Economy S V D to produce U sub H, Sigma sub H, and V sub H super T with rank r.

4. Branching Logic: The process splits into two paths based on rank r sub H.

* Path 1 (if r sub H equals r minus 1): Shows a partitioned matrix A sub 1 and B sub 1, leading to the Linear 1 model equations involving v sub H,1:r-1[k + 1] and a[k].

* Path 2 (else if r sub H equals r): Shows a single matrix A sub 2, leading to the Linear 2 model equations involving v sub H[k + 1] and a[k].

5. Bottom Section: Both paths converge to define parameters (A, C) and Q. These parameters, along with measurement noise g[k] and sensor input s[k], feed into a final box labeled Kalman filter. The filter contains equations for the state estimate v-hat, the output estimate a-hat, and the sensor estimate s-hat using the S U sub Y transformation.

2.3.1. First linear model based on time-delayed Hankel coefficients (linear 1)

We now define the vector $ {\boldsymbol{v}}_{H,1:r-1}\left({t}_j\right)={\left[{v}_{H,1}\left({t}_j\right),{v}_{H,2}\left({t}_j\right),\dots, {v}_{H,r-1}\left({t}_j\right)\right]}^{\top}\in {\mathrm{\mathbb{R}}}^{r-1} $ $ \left(j=1\dots p\right) $ extracted from the columns of matrix $ {\boldsymbol{V}}_H $ . We consider $ {\boldsymbol{v}}_{H,1:r-1} $ as the state variable and formulate the discrete in time, forced, linear dynamical system,

(2.12) $$ {\boldsymbol{v}}_{H,1:r-1}\left[k+1\right]={\boldsymbol{A}}_1{\boldsymbol{v}}_{H,1:r-1}\left[k\right]+{\boldsymbol{B}}_1{v}_{H,r}, $$

where $ {\boldsymbol{v}}_{H,1:r-1}\left[k\right]={\boldsymbol{v}}_{H,1:r-1}\left[{t}_k\right] $ , $ k=1\dots p-1 $ . We build a linear model using only the first ( $ r-1 $ ) elements of each column of $ {\boldsymbol{V}}_H $ , while the last element is used as forcing term.

The system matrices $ {\boldsymbol{A}}_1\in {\mathrm{\mathbb{R}}}^{\left(r-1\right)\times \left(r-1\right)} $ and $ {\boldsymbol{B}}_1\in {\mathrm{\mathbb{R}}}^{\left(r-1\right)} $ are computed by writing Eq. (2.12) explicitly for $ k=1 $ to $ p-1 $ ,

(2.13a) $$ {\boldsymbol{v}}_{H,1:r-1}\left[2\right]={\boldsymbol{A}}_1{\boldsymbol{v}}_{H,1:r-1}\left[1\right]+{\boldsymbol{B}}_1{v}_{H,r}\left[1\right], $$

(2.13b) $$ {\boldsymbol{v}}_{H,1:r-1}\left[3\right]={\boldsymbol{A}}_1{\boldsymbol{v}}_{H,1:r-1}\left[2\right]+{\boldsymbol{B}}_1{v}_{H,r}\left[2\right], $$

$$ \vdots $$

(2.13c) $$ {\boldsymbol{v}}_{H,1:r-1}\left[p\right]={\boldsymbol{A}}_1{\boldsymbol{v}}_{H,1:r-1}\left[p-1\right]+{\boldsymbol{B}}_1{v}_{H,r}\left[p-1\right]. $$

The above set can be written in matrix form as

(2.14) $$ {\boldsymbol{V}}_{1:r-1}^{\prime }={\boldsymbol{A}}_1{\boldsymbol{V}}_{1:r-1}+{\boldsymbol{B}}_1{\boldsymbol{V}}_r, $$

where $ {\boldsymbol{V}}_{1:r-1}^{\prime },{\boldsymbol{V}}_{1:r-1}\in {\mathrm{\mathbb{R}}}^{\left(r-1\right)\times \left(p-1\right)} $ and $ {\boldsymbol{V}}_r\in {\mathrm{\mathbb{R}}}^{p-1} $ take the following forms:

(2.15a) $$ {\boldsymbol{V}}_{1:r-1}^{\prime }={\left[{\boldsymbol{v}}_{H,1:r-1}\left[2\right]\hskip0.5em {\boldsymbol{v}}_{H,1:r-1}\left[3\right]\hskip0.5em \dots \hskip0.5em {\boldsymbol{v}}_{H,1:r-1}\left[p\right]\right]}^{\top }, $$

(2.15b) $$ {\boldsymbol{V}}_{1:r-1}={\left[\begin{array}{cccc}{\boldsymbol{v}}_{H,1:r-1}\left[1\right]& {\boldsymbol{v}}_{H,1:r-1}\left[2\right]& \dots & {\boldsymbol{v}}_{H,1:r-1}\left[p-1\right]\end{array}\right]}^{\top }, $$

(2.15c) $$ {\boldsymbol{V}}_r={\left[\begin{array}{cccc}{v}_{H,r}\left[1\right]& {v}_{H,r}\left[2\right]& \dots & {v}_{H,r}\left[p-1\right]\end{array}\right]}^{\top }. $$

Eq. (2.14) can be written in a more compact form as

(2.16a) $$ {\boldsymbol{V}}_{1:r-1}^{\prime }=\underset{=\boldsymbol{D}}{\underbrace{\left(\begin{array}{cc}{\boldsymbol{A}}_1& {\boldsymbol{B}}_1\end{array}\right)}}\underset{=\underline {\boldsymbol{V}}}{\underbrace{\left(\begin{array}{c}{\boldsymbol{V}}_{1:r-1}\\ {}{\boldsymbol{V}}_r\end{array}\right)}}, $$

where $ \boldsymbol{D}\in {\mathrm{\mathbb{R}}}^{\left(r-1\right)\times r} $ and $ \underline {\boldsymbol{V}}\in {\mathrm{\mathbb{R}}}^{r\times \left(p-1\right)} $ . In matrix $ \boldsymbol{D} $ , the first ( $ r-1 $ ) columns represent $ {\boldsymbol{A}}_1 $ and the last column represents $ {\boldsymbol{B}}_1 $ . So we have

(2.17) $$ {\boldsymbol{V}}_{1:r-1}^{\prime }=\boldsymbol{D}\underline {\boldsymbol{V}}\Rightarrow {\boldsymbol{V}}_{1:r-1}^{\prime }{\underline {\boldsymbol{V}}}^{\top }=\boldsymbol{D}\underline {\boldsymbol{V}}{\underline {\boldsymbol{V}}}^{\top}\Rightarrow \boldsymbol{D}={\boldsymbol{V}}_{1:r-1}^{\prime }{\underline {\boldsymbol{V}}}^{\top }{\left(\underline {\boldsymbol{V}}{\underline {\boldsymbol{V}}}^{\top}\right)}^{-1}. $$

Once $ {\boldsymbol{B}}_1 $ is known, the forcing $ {\boldsymbol{B}}_1{\boldsymbol{v}}_{H,r}\left[k\right] $ can be easily computed from the training dataset. The covariance matrix of the forcing $ {\boldsymbol{Q}}_1\in {\mathrm{\mathbb{R}}}^{\left(r-1\right)\times \left(r-1\right)} $ can be obtained from

(2.18) $$ {\boldsymbol{Q}}_1=\operatorname{cov}\left({\boldsymbol{B}}_1{v}_{H,r}\right)=\frac{1}{p-1}\sum \limits_{k=1}^{k=p-1}\left({\boldsymbol{B}}_1{v}_{H,r}\left[k\right]\right){\left({\boldsymbol{B}}_1{v}_{H,r}\left[k\right]\right)}^{\top }. $$

Vector $ {\boldsymbol{v}}_{H,1:r-1}\left[k\right] $ can be used to obtain the POD coefficients at the $ k $ -th instant from

(2.19) $$ \boldsymbol{a}\left[k\right]=\underset{={\boldsymbol{C}}_1}{\underbrace{{\boldsymbol{U}}_{H,1:r-1}\left({t}_1\right){\boldsymbol{\Sigma}}_{H,1:r-1}}}{\boldsymbol{v}}_{H,1:r-1}\left[k\right], $$

where $ {\boldsymbol{C}}_1\in {\mathrm{\mathbb{R}}}^{m_u\times \left(r-1\right)} $ and matrix $ {\boldsymbol{U}}_{H,1:r-1}\left({t}_1\right) $ consists of the $ \left(r-1\right) $ elements of the top row of $ {\boldsymbol{U}}_H $ , i.e.,

(2.20) $$ {\boldsymbol{U}}_{H,1:r-1}\left({t}_1\right)=\left[\begin{array}{cccc}{\boldsymbol{U}}_{H,1}^{\left(u,v,w\right)}\left({t}_1\right)& {\boldsymbol{U}}_{H,2}^{\left(u,v,w\right)}\left({t}_1\right)& \dots & {\boldsymbol{U}}_{H,r-1}^{\left(u,v,w\right)}\left({t}_1\right)\end{array}\right]. $$

2.3.2. Second linear model based on time-delayed Hankel coefficients (linear 2)

In the second linear model, the state variable is $ {\boldsymbol{v}}_H\left({t}_j\right)={\left[{v}_{H,1}\left({t}_j\right),{v}_{H,2}\left({t}_j\right),\dots, {v}_{H,r}\left({t}_j\right)\right]}^{\top}\in {\mathrm{\mathbb{R}}}^r $ $ \left(j=1\dots p\right) $ , i.e., we consider all $ r $ elements of each column of matrix $ {\boldsymbol{V}}_H $ . This model was first presented in Lu and Papadakis (Reference Lu and Papadakis2026a), but it is repeated below for completeness. A schematic overview of the algorithm is presented in Figure 1.

We seek a linear model of the form

(2.21) $$ {\boldsymbol{v}}_H\left[k+1\right]={\boldsymbol{A}}_2{\boldsymbol{v}}_H\left[k\right]+{\boldsymbol{w}}_2\left[k\right], $$

where $ {\boldsymbol{w}}_2\left[k\right] $ is random forcing. To obtain $ {\boldsymbol{A}}_2\in {\mathrm{\mathbb{R}}}^{r\times r} $ , we write Eq. (2.21) explicitly for $ k=1 $ to $ p-1 $ :

(2.22a) $$ {\boldsymbol{v}}_H\left[2\right]\approx {\boldsymbol{A}}_2{\boldsymbol{v}}_H\left[1\right], $$

(2.22b) $$ {\boldsymbol{v}}_H\left[3\right]\approx {\boldsymbol{A}}_2{\boldsymbol{v}}_H\left[2\right], $$

$$ \vdots $$

(2.22c) $$ {\boldsymbol{v}}_H\left[p\right]\approx {\boldsymbol{A}}_2{\boldsymbol{v}}_H\left[p-1\right]. $$

In matrix form, this becomes

(2.23) $$ {\boldsymbol{V}}^{\prime}\approx {\boldsymbol{A}}_2\boldsymbol{V}, $$

where $ {\boldsymbol{V}}^{\prime },\boldsymbol{V}\in {\mathrm{\mathbb{R}}}^{r\times \left(p-1\right)} $ .

(2.24a) $$ {\boldsymbol{V}}^{\prime }={\left[\begin{array}{cccc}{\boldsymbol{v}}_H\left[2\right]& {\boldsymbol{v}}_H\left[3\right]& \dots & {\boldsymbol{v}}_H\left[p\right]\end{array}\right]}^{\top }, $$

(2.24b) $$ \boldsymbol{V}={\left[\begin{array}{cccc}{\boldsymbol{v}}_H\left[1\right]& {\boldsymbol{v}}_H\left[2\right]& \dots & {\boldsymbol{v}}_H\left[p-1\right]\end{array}\right]}^{\top }. $$

Then, system matrix $ {\boldsymbol{A}}_2 $ can be calculated by

(2.25) $$ {\boldsymbol{V}}^{\prime }={\boldsymbol{A}}_2\boldsymbol{V}\Rightarrow {\boldsymbol{V}}^{\prime }{\boldsymbol{V}}^{\top }={\boldsymbol{A}}_2{\boldsymbol{V}\boldsymbol{V}}^{\top}\Rightarrow {\boldsymbol{A}}_2={\boldsymbol{V}}^{\prime }{\boldsymbol{V}}^{\top }{\left({\boldsymbol{V}\boldsymbol{V}}^{\top}\right)}^{-1}. $$

Once $ {\boldsymbol{A}}_2 $ is known, the forcing $ {\boldsymbol{w}}_2\left[k\right] $ can be easily calculated from the training dataset, $ {\boldsymbol{w}}_2\left[k\right]={\boldsymbol{v}}_H\left[k+1\right]-{\boldsymbol{A}}_2{\boldsymbol{v}}_H\left[k\right] $ . As before, the covariance matrix $ {\boldsymbol{Q}}_2\in {\mathrm{\mathbb{R}}}^{r\times r} $ of $ {\boldsymbol{w}}_2 $ can be obtained from,

(2.26) $$ {\boldsymbol{Q}}_2=\unicode{x1D53C}\left({\boldsymbol{w}}_2{\boldsymbol{w}}_2^{\top}\right)=\frac{1}{p-1}\sum \limits_{k=1}^{k=p-1}{\boldsymbol{w}}_2\left[k\right]{\boldsymbol{w}}_2^{\top}\left[k\right]. $$

From vector $ {\boldsymbol{v}}_H\left[k\right] $ , we can get

(2.27) $$ \boldsymbol{a}\left[k\right]=\underset{={\boldsymbol{C}}_2}{\underbrace{{\boldsymbol{U}}_H\left({t}_1\right){\boldsymbol{\Sigma}}_H}}{\boldsymbol{v}}_H\left[k\right], $$

where $ {\boldsymbol{C}}_2\in {\mathrm{\mathbb{R}}}^{m_u\times r} $ and $ {\boldsymbol{U}}_H\left({t}_1\right) $ is the top row of $ {\boldsymbol{U}}_H $ .

2.3.3. N4sid algorithm

A schematic overview of this algorithm is displayed in Figure 2. In n4sid, we seek a linear dynamical system of the form

(2.28) $$ {\displaystyle \begin{array}{ll}\boldsymbol{x}\left[k+1\right]& =\underset{\mathrm{System}\ \mathrm{matrix}}{\underbrace{\boldsymbol{A}}}\boldsymbol{x}\left[k\right]+\underset{\mathrm{noise}}{\underbrace{\boldsymbol{w}\left[k\right]}}\\ {}\underset{\mathrm{Output}}{\underbrace{\boldsymbol{a}\left[k\right]}}& =\underset{\mathrm{Output}\ \mathrm{matrix}}{\underbrace{\boldsymbol{C}}}\boldsymbol{x}\left[k\right]+\underset{\mathrm{noise}}{\underbrace{\boldsymbol{v}\left[k\right]}},\end{array}} $$

where $ \boldsymbol{a}\left[k\right]=\boldsymbol{a}\left({t}_K\right) $ , $ \boldsymbol{x}\in {\mathrm{\mathbb{R}}}^n $ is the internal state vector of size $ n $ , $ \boldsymbol{A}\in {\mathrm{\mathbb{R}}}^{n\times n} $ and $ \boldsymbol{C}\in {\mathrm{\mathbb{R}}}^{m_u\times n} $ , $ \boldsymbol{w}\left[k\right]\in {\mathrm{\mathbb{R}}}^n $ and $ \boldsymbol{v}\left[k\right]\in {\mathrm{\mathbb{R}}}^{m_u} $ are Gaussian white noise sequences. Note that we have introduced the internal state $ \boldsymbol{x}\left[k\right] $ with size $ n $ , which is independent of the total number of retained POD modes, $ {m}_u $ . This gives additional flexibility, as it allows one to adjust $ n $ and thus the size of the matrices $ \boldsymbol{A} $ and $ \boldsymbol{C} $ , so that (2.28) returns a vector sequence $ \boldsymbol{a}\left[k\right] $ that optimally fits the true values.

Figure 2.

Schematic overview of the algorithm based on the system identification method n4sid.

Figure 2. Long description

The flowchart is organized into three main horizontal stages.

Top row: Flow field data is shown as a stack of three heatmaps labeled k equals 1, k equals 2, and k equals K. An arrow points to an Economy S V D section where a large blue matrix V super 1/2 Y is approximately decomposed into three smaller matrices: U sub Y, Sigma sub Y (a diagonal matrix in a red dashed box), and V sub Y super T.

Middle row: An arrow from V sub Y super T leads to a set of four time-series plots labeled True output: a[k] for m equals 1 through 4. These plots feed into an n 4 s i d block, which outputs a Plant state-space model: x[k plus 1] equals A x[k] plus w[k] and a[k] equals C x[k] plus v[k]. Below the true output, four heatmaps labeled U sub Y, 1 and U sub Y, 2 for variables u prime and v prime are multiplied by the true output to produce S U sub Y.

Bottom row: Sensor measurements are shown as two line graphs with red and blue waves. These provide the sensor input s[k]. The plant model parameters (A, C) and Q, along with measurement noise g[k] and R, feed into a Kalman filter block. The filter contains the equations: x-hat[k plus 1] equals A x-hat[k] plus L (s[k] minus s-hat[k]) and a-hat[k] equals C x-hat[k].

We use the system identification algorithm n4sid to obtain $ \boldsymbol{A} $ and $ \boldsymbol{C} $ , as well as the covariance matrices of sequences $ \boldsymbol{w}\left[k\right] $ and $ \boldsymbol{v}\left[k\right] $ . The algorithm also assembles the time-delayed Hankel matrix (2.8) and considers the top $ q/2 $ rows as past outputs and the bottom $ q/2 $ as future outputs. The next step is the orthogonal projection of the row space of the future outputs into the row space of the past outputs and the SVD of the resulting matrix. In many cases, there is only a small number of large singular values and the rest can be discarded; this yields an optimal order $ n $ . Alternatively, the order $ n $ can be prescribed. From the left and right singular vectors, the states $ \boldsymbol{x}\left[k\right] $ are obtained, and then, the matrices $ \boldsymbol{A},\boldsymbol{C} $ are computed using least squares. Note that $ \boldsymbol{A} $ and $ \boldsymbol{C} $ are unique up to a similarity transformation. More details can be found in Van Overschee and De Moor (Reference Van Overschee and De Moor1994) and the book of Van Overschee and De Moor (Reference Van Overschee and De Moor1996), in particular chapter 3 that deals with stochastic identification (i.e., when only output data are available, as in our case). The algorithm is implemented in MATLAB and can be executed with the command n4sid.

In conclusion, n4sid also starts with the same Hankel matrix (2.8) as the previous two algorithms, but then follows a different path. The states $ \boldsymbol{x}\left[k\right] $ are obtained directly from the output vectors, and then the system matrices are computed. On the other hand, for linear 1 and linear 2, the dominant modes of (2.8) are first determined, and then, a system is constructed from the evolution of the modal coefficients.

2.4. Closure of the systems using sensor measurements

The linear systems (2.12), (2.21), or (2.28) cannot be used directly because the forcing terms are unknown. However, if sensor measurements are available, they can be used to close the systems. This is achieved as explained below.

Let us assume that we have a set of $ l $ velocity measurements $ \boldsymbol{s}\left[k\right]\in {\mathrm{\mathbb{R}}}^l $ at a number of sensor points. At each sensor, one or more velocity components are recorded. We can express $ \boldsymbol{s}\left[k\right] $ in terms of $ \boldsymbol{a}\left[k\right] $ as follows

(2.29) $$ \boldsymbol{s}\left[k\right]={\boldsymbol{SU}}_Y\boldsymbol{a}\left[k\right]+\boldsymbol{g}\left[k\right], $$

where matrix $ \boldsymbol{S} $ selects the rows of the POD mode matrix $ {\boldsymbol{U}}_Y $ (see (2.7)) corresponding to the sensor location and the velocity component being measured (it is $ 0 $ everywhere except for those points and components where the corresponding element is equal to $ 1 $ ). Vector $ \boldsymbol{g}\in {\mathrm{\mathbb{R}}}^l $ includes measurement errors as well as errors due to POD mode truncation (because only $ {m}_u $ modes are retained). Here, we assume that $ \boldsymbol{g}\left[k\right] $ arises only due to the POD truncation; thus, the elements of $ \boldsymbol{g}\left[k\right] $ can be obtained from the training dataset, $ \boldsymbol{g}\left[k\right]=\boldsymbol{s}\left[k\right]-{\boldsymbol{SU}}_Y\boldsymbol{a}\left[k\right] $ . However, this is not a restriction; sensor noise can be easily included. The covariance $ \boldsymbol{R}\in {\mathrm{\mathbb{R}}}^{l\times l} $ can be easily calculated from

(2.30) $$ \boldsymbol{R}=\unicode{x1D53C}\left({\boldsymbol{g}\boldsymbol{g}}^{\top}\right)=\frac{1}{p-1}\sum \limits_{k=1}^{k=p-1}\boldsymbol{g}\left[k\right]{\boldsymbol{g}}^{\top}\left[k\right]. $$

We are now ready to design a Kalman filter to estimate $ {\boldsymbol{v}}_{H,1:{r}_H}\left[k\right] $ from the sparse measurements $ \boldsymbol{s}\left[k\right] $ , where $ {r}_H $ is either $ \left(r-1\right) $ or $ r $ depending on whether the forcing uses the $ r $ -th element (linear 1) or it is random (linear 2), respectively. In either case, the filter takes the form

(2.31a) $$ {\hat{\boldsymbol{v}}}_{H,1:{r}_H}\left[k+1\right]=\boldsymbol{A}{\hat{\boldsymbol{v}}}_{H,1:{r}_H}\left[k\right]+\mathrm{\mathcal{L}}\left(\boldsymbol{s}\left[k\right]-\hat{\boldsymbol{s}}\left[k\right]\right), $$

(2.31b) $$ \hat{\boldsymbol{a}}\left[k\right]=\boldsymbol{C}{\hat{\boldsymbol{v}}}_{H,1:{r}_H}\left[k\right], $$

(2.31c) $$ \hat{\boldsymbol{s}}\left[k\right]={\boldsymbol{SU}}_Y\hat{\boldsymbol{a}}\left[k\right]={\boldsymbol{SU}}_Y\boldsymbol{C}{\hat{\boldsymbol{v}}}_{H,1:{r}_H}\left[k\right], $$

The hat $ \hat{\left(\cdot \right)} $ denotes an estimated quantity and $ \mathrm{\mathcal{L}}\in $ $ {\mathrm{\mathbb{R}}}^{r_H\times l} $ is the Kalman filter gain. The latter is obtained from the solution of the following Riccati equation:

(2.32a) $$ \boldsymbol{P}={\boldsymbol{APA}}^T-\boldsymbol{AP}{\left({\boldsymbol{SU}}_Y\boldsymbol{C}\right)}^T{\left({\boldsymbol{SU}}_Y\boldsymbol{C}\boldsymbol{P}{\left({\boldsymbol{SU}}_Y\boldsymbol{C}\right)}^T+\boldsymbol{R}\right)}^{-1}{\boldsymbol{SU}}_Y{\boldsymbol{CPA}}^T+\boldsymbol{Q}, $$

(2.32b) $$ \mathrm{\mathcal{L}}=\boldsymbol{AP}{\left({\boldsymbol{SU}}_Y\boldsymbol{C}\right)}^T{\left({\boldsymbol{SU}}_Y\boldsymbol{C}\boldsymbol{P}{\left({\boldsymbol{SU}}_Y\boldsymbol{C}\right)}^T+\boldsymbol{R}\right)}^{-1}. $$

We use the same idea to close the n4sid system (2.28).

6.1. Assessment of reconstruction quality of different linear models

We now compare the reconstruction quality of the three estimators based on models linear 1, linear 2, and n4sid. The former two have different forcing terms, as explained in Section 2. Due to this difference, linear 1 has order ( $ r-1 $ ) and linear 2 has order $ r $ , while n4sid has order $ n $ .

The reconstruction quality is quantified with the metric $ \mathrm{FIT}\left[\%\right] $ defined as

(6.1) $$ \mathrm{FIT}\left[\%\right]=100\left(1-\frac{\sum_{i=1}^{m_u}\overline{{\left({a}_i\left[k\right]-{\hat{a}}_i\left[k\right]\right)}^2}}{\sum_{i=1}^{m_u}\overline{{\left({a}_i\left[k\right]\right)}^2}}\right), $$

where $ {m}_u $ is the number of retained velocity POD modes. $ \mathrm{FIT}\left[\%\right] $ measures how close the estimated velocity POD coefficients track the true coefficients. For perfect matching, it is equal to 100%, but can become negative for large deviations between the true and estimated values.

We start with time-delayed embedding $ q=1 $ and use $ {m}_u=7 $ velocity modes and $ {m}_c=9 $ scalar modes (that capture 97.46% of the energy of the fluctuating velocity field and 97.43% of the scalar variance). For linear 1 and linear 2, the maximum order $ {r}_{\mathrm{max}} $ is set by the number of modes used and the time delay; it is equal to $ {r}_{\mathrm{max}}={m}_u\times q $ when reconstructing the flow using velocity measurements or $ {r}_{\mathrm{max}}=\left({m}_u+{m}_c\right)\times q $ when reconstructing from scalar measurements. We set $ r $ to the maximum value, i.e., $ r={m}_u=7 $ or $ r={m}_u+{m}_c=16 $ (because we consider $ q=1 $ ). On the other hand, the order $ n $ of n4sid is set independently from $ {m}_u $ and $ {m}_c $ because it employs an internal state (equation (2.28)).

In Figure 4, the reconstruction quality of the three linear models is compared. For a meaningful comparison, we consider similar values of $ r $ and $ n $ . In panel ( $ a $ ), the n4sid model with orders $ n=\mathrm{7,8,9,10} $ is compared with the two linear models with $ r={m}_u=7 $ (reconstruction from velocity sensors). The performance of n4sid is poor for $ n=7 $ , but improves significantly for a larger $ n $ , and for $ n=10 $ , the results become almost identical to linear 2. Notice that even with one sensor point, all models can provide FIT[ $ \% $ ] values more than 90%. Panel ( $ b $ ) refers to reconstruction from scalar measurements, and the order of n4sid is increased to $ n=\mathrm{20,25,30,35} $ . Higher $ n $ improves the reconstruction quality, but it is not until $ n=35 $ that n4sid yields a reconstruction quality similar to that of the two other models. Moreover, the linear 2 model performs significantly better than the linear 1 model. For a single scalar sensor, the FIT[ $ \% $ ] values are between 80 and 90%, slightly smaller than velocity sensors. This is probably because a scalar sensor provides only a single value, while a velocity sensor records both $ u $ and $ v $ velocity components. The conclusion from this figure is that all models deliver similar results, but n4sid requires higher order. It is also much more time-consuming. For $ n=35 $ , n4sid takes more than $ 30 $ minutes on a desktop with Intel Core i7 9700 3.02666 MHz and 32GB DIMM memory, while for $ r=16 $ , linear 1 and linear 2 take only less than $ 1 $ minute and give comparable results. Linear 1 and linear 2 are therefore more robust and computationally efficient.

Figure 4.

Reconstruction quality FIT[ $ \% $ ] of different linear models against the number of ( $ a $ ) velocity sensors and ( $ b $ ) scalar sensors.

Figure 4. Long description

Two side-by-side line graphs, labeled a and b, share a common y-axis of F I T percentage ranging from 40 to 100 and an x-axis of Number of sensors, n sub p, ranging from 1 to 16.

Panel a, velocity sensors:

* Most models, including n 4 s i d with n equals 8, 9, 10, and linear 1 and 2 with r equals 7, show high, stable performance between 95 and 100 percent F I T.

* The n 4 s i d, n equals 7 model is an outlier, starting at 40 percent and rising to a plateau around 65 percent.

Panel b, scalar sensors:

* Performance is more varied. The n 4 s i d, n equals 35 and linear 2, r equals 16 models maintain the highest accuracy, plateauing near 98 percent.

* Linear 1, r equals 16 and n 4 s i d, n equals 30 follow, plateauing between 90 and 95 percent.

* The n 4 s i d models with n equals 20 and 25 show significant volatility for low sensor counts, eventually stabilizing between 80 and 88 percent.

Legends in the bottom-right of each panel identify the models by color and marker shape: red squares, black circles, and blue diamonds, with solid lines for lower n or r values and dashed lines for higher values.

We now compare the performance of the three models from scalar measurements when the same order is used. We set the model order of n4sid to $ n=\mathrm{20,25,30} $ . For linear 1 and linear 2, we need to increase the time delay to $ q=2 $ to match $ n $ ; this results in $ {r}_{max}=\left({m}_u+{m}_c\right)\times q=32 $ , and when performing the SVD of the Hankel matrix in (3.6), we truncate $ r $ so that $ r=n $ . The results from this comparison are shown in Figure 5. It can be seen that linear 1 and linear 2 models give consistently better reconstruction than n4sid model for all values of $ n $ . As $ r $ increases, the performance of linear 1 and linear 2 improves, and the discrepancies between them decrease. The linear 2 model slightly outperforms the linear 1 model. For $ {n}_p<5 $ , n4sid gives a poor reconstruction quality, but as $ {n}_p $ increases, the quality gradually improves and reaches a plateau. Even with the largest $ {n}_p $ , the performance of the n4sid model is still lower than that of the other two models. Again, linear 1 and linear 2 models significantly outperform the n4sid model and achieve a high level of reconstruction quality with only one or two sensors.

Figure 5.

FIT[%] against the number of scalar sensors, $ {n}_p $ , for different model orders $ n=r $ , using the Kalman filter based on the n4sid model (□), the linear 1 model (□), and the linear 2 model (□).

Figure 5. Long description

A three-panel line graph labeled a, b, and c. All panels share a common Y-axis labeled F I T in percent ranging from 30 to 100 and an X-axis labeled Number of sensors, n sub p, ranging from 1 to 16.

* Panel a, n equals r equals 20. The blue linear 2 line starts at 87 percent and stabilizes near 98 percent. The black linear 1 line starts at 88 percent and stabilizes near 96 percent. The red n 4 s i d line starts at 31 percent, rises sharply to 78 percent at 4 sensors, and fluctuates between 80 and 84 percent.

* Panel b, n equals r equals 25. The blue linear 2 line starts at 90 percent and reaches 99 percent. The black linear 1 line starts at 87 percent and reaches 96 percent. The red n 4 s i d line starts at 50 percent, drops to 38 percent at 3 sensors, then climbs to a plateau around 90 percent.

* Panel c, n equals r equals 30. The blue linear 2 line starts at 91 percent and reaches 99 percent. The black linear 1 line starts at 90 percent and reaches 98 percent. The red n 4 s i d line starts at 49 percent, climbs to 79 percent at 5 sensors, and reaches a peak of 86 percent at 9 sensors before stabilizing around 85 percent.

In all panels, linear 2 consistently performs best, followed closely by linear 1, while n 4 s i d shows the lowest accuracy and highest sensitivity to the number of sensors.

A possible explanation for the better performance of these two algorithms is the following. By construction, they extract the most dominant modes of the Hankel matrix and derive a linear model that tracks the evolution of the modal coefficients. On the other hand, n4sid first derives a mapping from past to future outputs without considering the underlying structure of the outputs (recall that the outputs are the POD coefficients). Perhaps a rescaling of the outputs, e.g., taking into account the energy of each POD mode, would improve the performance.

The better reconstruction quality of linear 2 compared to linear 1 can be explained by inspecting the probability density function (PDF) of the last element $ {v}_{H,r} $ that acts as a forcing in the linear 1 model (see (2.12)). This is shown in Figure 6; as can be seen, the PDF is very close to Gaussian. For the three-equation Lorenz system, $ {v}_r\left[k\right] $ is intermittent (is activated only prior to lobe switching), and therefore, its PDF is far from Gaussian (see SL Brunton et al., Reference Brunton, Brunton, Proctor, Kaiser and Kutz2017). In the present case, the Hankel matrix consists of a finite number of $ {m}_u $ (or $ {m}_u+{m}_c $ ) POD coefficients. The set of governing equations for this truncated set can be found in Guzmán Iñigo et al. (Reference Guzmán Iñigo, Sodar and Papadakis2019). In this set, the forcing term consists of a linear part (due to truncation) and a nonlinear part (due to a quadratic interaction); refer to equations 9(a) and 9(b) in the aforementioned reference. It is well known that high-order POD coefficients exhibit Gaussian behavior (Berkooz et al., Reference Berkooz, Holmes and Lumley1993), so if the linear forcing term dominates over the nonlinear, this may explain the Gaussian distribution of $ {v}_r $ . For the time-delayed linear 1 model, $ {v}_r $ is treated as a forcing term, but since it is Gaussian, it does not offer much additional information (as opposed, e.g., to the Lorenz system). Therefore, the linear 2 model that has a slightly higher order than linear 1 ( $ r $ instead of $ r-1 $ ) provides slightly better results.

Figure 6.

Probability density function (PDF) of the forcing term $ {v}_{H,16} $ (red), plotted against an appropriately scaled normal distribution (black).

Figure 7 shows the true and reconstructed vorticity fields at two time instants using only two scalar sensors. The main features, such as the large-scale vortices, are well reproduced in both instances by all models. However, the n4sid model (second row) gives spurious results in the near wake, but linear 1 and linear 2 models are much better at reconstructing the vortices close to the prism. n4sid results also show regions of unexpected positive or negative vorticity far downstream; again, this indicates the lack of robustness.

Figure 7.

Contour plots of instantaneous vorticity: (left column) at $ t=3.76 $ , (right column) at $ t=11.36 $ . ( $ a $ - $ b $ , top row) DNS data, reconstruction using the Kalman filter based on ( $ c $ - $ d $ , second row) the n4sid model with $ n=30 $ , ( $ e $ - $ f $ , third row) linear 1 model, and ( $ g $ - $ h $ , bottom row) linear 2 model with $ r=30 $ . All results are obtained with $ {n}_p=2 $ . Square (■) markers indicate the location of the $ {c}^{\prime } $ sensors.

Figure 7. Long description

A multi-panel display of eight contour plots arranged in a 4 by 2 grid. The x-axis for all plots ranges from negative 2 to 10, labeled x all over h. The y-axis ranges from 0 to 4, labeled y all over h. A gray rectangular block is positioned at the origin on the bottom left of each plot.

* The left column represents time t = 3.76 and the right column represents time t = 11.36.

* Row 1 (panels a and b) shows D N S data. Panel a shows a concentrated blue vortex at x = 9 and a smaller red vortex at x = 7. Panel b shows a large blue vortex centered at x = 4.

* Row 2 (panels c and d) shows the n 4 s i d model with n = 30. Two black square sensors labeled 1 and 2 are located at x = 2.5. The contours are more fragmented compared to the D N S data.

* Row 3 (panels e and f) shows the linear 1 model. It captures the primary vortex structures but with smoother, less detailed gradients than the D N S.

* Row 4 (panels g and h) shows the linear 2 model with r = 30. The results are visually similar to the linear 1 model, showing the evolution of the wake behind the block.

Each row includes a vertical color scale on the right ranging from negative 6 (dark blue) to 0 (white) to 6 (dark red), representing the vorticity magnitude.

The above observations can be explained by inspecting the POD coefficients of the seven most dominant velocity modes (Figure 8). The estimated time coefficients of linear 1 and linear 2 models (middle and right columns) start from zero initial conditions and quickly catch up with the true time coefficients. On the other hand, the n4sid model (left column) produces a phase shift between the true and estimated coefficients of the first and second modes. This explains the slight spatial shift in the location of the main vortex in Figure 7. The n4sid model also overestimates the time coefficients of higher modes, and this may explain the spurious features in the second row of Figure 7. The linear 1 and linear 2 models provide much better reconstruction quality.

Figure 8.

POD time coefficients $ {a}_1(t) $ to $ {a}_7(t) $ for velocity modes on the validation dataset: reconstruction using the Kalman filter based on the (left column $ a $ ) n4sid model with $ n=30 $ , (middle column $ b $ ) linear 1 model, (right column $ c $ ) linear 2 model with $ r=30 $ . All results are obtained with $ {n}_p=2 $ . DNS: blue solid line. Reconstruction from the model with zero initial conditions: red dashed line.

Figure 8. Long description

The multi-panel line graph consists of 21 individual plots arranged in three columns labeled a, b, and c, and seven rows labeled a sub 1 (t) through a sub 7 (t). The x-axis for all plots is Time, ranging from 0 to approximately 95. The y-axis represents the coefficient value, with scales decreasing in magnitude from the top row to the bottom row.

In every plot, a blue solid line represents the D N S ground truth, showing periodic oscillations. A red dashed line represents the model reconstruction.

* Column a (n 4 s i d model): The red dashed line shows significant phase and amplitude discrepancies compared to the blue line, particularly in the lower rows (a sub 5 through a sub 7), where the reconstruction appears noisy and poorly aligned.

* Column b (linear 1 model): The red dashed line tracks the blue solid line much more closely than in column a. There is a slight initial transient period at the start of the time series, after which the lines synchronize well across all seven modes.

* Column c (linear 2 model): This column shows the highest level of agreement. The red dashed line overlaps almost perfectly with the blue solid line from the beginning of the time series through to the end, maintaining accuracy even in the high-frequency oscillations of the lower modes.

We set $ n=r=30 $ and $ {n}_p=2 $ (scalar sensors) and compare the flow statistics (Reynolds stresses and kinetic energy of the fluctuating field) against the DNS results in Figure 9. The three models give similar results for $ \left\langle {u}^{\prime }{u}^{\prime}\right\rangle $ (left column), while n4sid model predicts spurious regions of high $ \left\langle {v}^{\prime }{v}^{\prime}\right\rangle $ downstream (middle column). In contrast, linear 1 and linear 2 models provide better reconstruction of $ \left\langle {v}^{\prime }{v}^{\prime}\right\rangle $ . The linear 1 model (third row) underestimates the $ \left\langle {v}^{\prime }{v}^{\prime}\right\rangle $ behind the prism, while the linear 2 model does better, confirming its superiority. Similar observations can be made for the reconstruction quality of the fluctuating kinetic energy (right column).

Figure 9.

Contour plots of flow statistics: (left column) $ \left\langle {u}^{\prime }{u}^{\prime}\right\rangle $ , (middle column) $ \left\langle {v}^{\prime }{v}^{\prime}\right\rangle $ , (right column) kinetic energy of the fluctuating field $ k $ . ( $ a $ - $ c $ , top row) DNS data, reconstruction using the Kalman filter based on ( $ d $ - $ f $ , second row) the n4sid model with $ n=30 $ , ( $ g $ - $ i $ , third row) linear 1 model and ( $ j $ - $ l $ , bottom row) linear 2 model with $ r=30 $ . All results are obtained with $ {n}_p=2 $ scalar sensors, indicated by square (■) markers.

Figure 9. Long description

The figure consists of 12 panels labeled a through l. Each panel shows a flow field with a gray square obstacle at the origin. The x-axis represents x over h from negative 2 to 10, and the y-axis represents y over h from 0 to 4.

* Columns: The left column (a, d, g, j) shows the expected value of u prime u prime. The middle column (b, e, h, k) shows the expected value of v prime v prime. The right column (c, f, i, l) shows the kinetic energy k.

* Rows: The top row (a-c) displays D N S data. The second row (d-f) shows the n 4 s i d model with n equals 30. The third row (g-i) shows the linear 1 model. The bottom row (j-l) shows the linear 2 model with r equals 30.

* Features: In rows 2 through 4, two black square sensors labeled 1 and 2 are positioned downstream of the obstacle at x over h equals 2.5.

* Color Scales: Each column has a color bar at the top. The left and middle columns range from 0 (dark blue) to 0.9 (dark red). The right column ranges from 0 to 0.54.

* Trends: The D N S data shows high-intensity turbulent structures shedding from the obstacle. The n 4 s i d model (second row) captures the spatial distribution and intensity most accurately compared to the D N S. The linear 1 and linear 2 models (rows 3 and 4) show smoother, less intense reconstructions of the fluctuating fields.

6.2. Forecasting the velocity field from scalar measurements

As mentioned in Section 4, we can forecast the evolution of the future velocity field from measurements at the current time instant. To this end, we increase the time-delay embedding parameter to $ q=\mathrm{400,800,1200} $ and $ 1600 $ with corresponding time delays $ q\times \Delta t=\mathrm{16,32,48} $ and $ 64 $ time units. The smallest time delay $ 16 $ is close to the period $ 15.84 $ of the first two velocity POD modes. The training dataset is $ {K}_{\mathrm{train}}\times \Delta t=95 $ . We set $ r=20 $ and $ {n}_p=1 $ to test the forecasting performance with only one scalar sensor.

The six most dominant left singular vectors $ {\boldsymbol{U}}_{H,i}^{\left(u,v\right)}\in {\mathrm{\mathbb{R}}}^{m_u} $ and $ {\boldsymbol{U}}_{H,i}^{(c)}\in {\mathrm{\mathbb{R}}}^{m_c} $ (see (3.7)) for $ {m}_u=7 $ , $ {m}_c=9 $ and $ q=1200 $ , $ q\times \Delta t=48 $ are shown in Figure 10. The horizontal axis is the time-delay $ j\times \Delta t $ , and the vertical axis the POD mode index $ k\left(=1\dots {m}_u\right) $ or $ l\left(=1\dots {m}_c\right) $ . It can be seen that the periodic behavior is imprinted in the structure of the singular modes of the time-delayed Hankel matrix. The two most dominant singular vectors fit three shedding periods (as expected for the considered $ q\times \Delta t $ ), and the structures are time-shifted to account for the vortex propagation once it is released in the wake. Higher-order singular modes capture higher harmonics. In Figure 11, we plot the future evolution of the first coefficient of the velocity and scalar modes, $ {a}_1(t) $ and $ {b}_1(t) $ , respectively, for different forecasting windows $ q\times \Delta t $ . These results are obtained from streaming measurements from a single scalar sensor; the last record is at $ t=0 $ , and only the predicted future evolution for $ t>0 $ is shown. As can be seen, $ {a}_1(t) $ and $ {b}_1(t) $ are reproduced excellently. Both coefficients are periodic with the same period, but of course, the time signals are different (notice for example that the amplitude of $ {a}_1(t) $ is five times larger than $ {b}_1(t) $ ). The algorithm learns the relative strengths (and phase difference) between the velocity and scalar fields, and once streaming data (even from a single sensor) become available, it can provide very reliable predictions of the future evolution. We have looked at the forecasting of the other modes; the quality slightly degrades for higher modes, but the $ \mathrm{FIT}\left[\%\right] $ metric is still above $ 60\% $ for the seventh mode (results not shown for brevity).

Figure 10.

Contours of the six most dominant left singular vectors of the time-delayed Hankel matrix in the time-delay/mode index plane for $ q=1200,q\times \Delta t=48 $ : (left column) $ {\mathbf{U}}_{H,i}^{\left(u,v\right)} $ of velocity POD modes (right column) $ {\mathbf{U}}_{H,i}^{(c)} $ of scalar POD modes.

Figure 10. Long description

The multi-panel figure consists of 12 horizontal contour plots arranged in two columns of six.

Common Axes and Scales:

- The horizontal x-axis for all plots is labeled j times Delta t, ranging from 0 to 48.

- The left column (velocity modes) has a y-axis labeled k, ranging from 1 to 7.

- The right column (scalar modes) has a y-axis labeled l, ranging from 1 to 9.

- Each plot includes a color bar on the right, where red indicates positive values and blue indicates negative values.

- Each plot is labeled with an index i from 1 to 6 in the top-left corner.

Left Column (Velocity P O D Modes):

- i equals 1 and 2: Show broad, elongated structures concentrated between k equals 1 and 3, with a slight tilt.

- i equals 3 and 4: Display a more frequent, periodic wave-like pattern of alternating red and blue diagonal bands extending up to k equals 5.

- i equals 5 and 6: Feature high-frequency, vertical zig-zag patterns with smaller spatial scales.

Right Column (Scalar P O D Modes):

- i equals 1 and 2: Show very smooth, low-intensity structures near the bottom of the l-axis.

- i equals 3 and 4: Exhibit distinct, alternating red and blue diagonal wave packets similar to the velocity modes but spanning l equals 1 to 5.

- i equals 5 and 6: Show complex, high-frequency oscillations that appear more fragmented than the velocity counterparts, spanning up to l equals 7.

Figure 11.

Forecasting of the POD temporal coefficients (left column) $ {a}_1(t) $ and (right column) $ {b}_1(t) $ : prediction using Kalman filter based on the time-delayed Hankel matrix of ( $ a $ – $ b $ ) $ q\times \Delta t=16 $ , ( $ c $ – $ d $ ) $ q\times \Delta t=32 $ , ( $ e $ – $ f $ ) $ q\times \Delta t=48 $ and ( $ g $ – $ h $ ) $ q\times \Delta t=64 $ . DNS: blue solid line. Model prediction: red dashed line.

Figure 11. Long description

The figure consists of eight panels arranged in a four-by-two grid. The left column (panels a, c, e, g) displays the temporal coefficient a sub 1 (t) on the y-axis, ranging from -0.5 to 0.5. The right column (panels b, d, f, h) displays b sub 1 (t) on the y-axis, ranging from -0.1 to 0.1. The x-axis represents Time for all panels.

* Row 1 (panels a and b): Time scale 0 to 16. Both graphs show a single wave cycle starting high, dipping to a trough around time 8, and rising again.

* Row 2 (panels c and d): Time scale 0 to 32. The graphs show two complete oscillatory cycles.

* Row 3 (panels e and f): Time scale 0 to 48. The graphs show three complete oscillatory cycles.

* Row 4 (panels g and h): Time scale 0 to 64. The graphs show four complete oscillatory cycles.

In every panel, a blue solid line representing D N S data is almost perfectly overlaid by a red dashed line representing the model prediction. The red dashed line tracks the peaks, troughs, and slopes of the blue line with high precision across all time intervals, indicating high model accuracy.

6.3. Flow field reconstruction from scalar measurements at the off-design condition Re = 800

The robustness of the three models designed at $ \mathit{\operatorname{Re}}=1000 $ is now assessed at the off-design condition of $ \mathit{\operatorname{Re}}=800 $ . The time coefficients at the new $ \mathit{\operatorname{Re}} $ are obtained by projecting the velocity and scalar snapshot data to the POD modes of $ \mathit{\operatorname{Re}}=1000 $ . Therefore, the models are trained at one operating condition and then applied to another; this makes them truly predictive.

The $ \mathrm{FIT}\left[\%\right] $ metric against the number of (scalar) sensors for three different model orders is shown in Figure 12. This is equivalent to Figure 5 for the reference case of $ \mathit{\operatorname{Re}}=1000 $ . As expected, the reconstruction quality at the off-design condition is not as good as for the baseline condition. Nevertheless, again linear 1 and linear 2 models offer better and more robust performance compared to the n4sid model with very small number of sensors.

Figure 12.

$ FIT\left[\%\right] $ against the number of scalar sensors, $ {n}_p $ , for different model orders $ n=r $ , using Kalman filter based on the n4sid algorithm (□), the linear 1 model (□) and the linear 2 model (□) at the off-design condition of $ \mathit{\operatorname{Re}}=800 $ .

Figure 12. Long description

The multi-panel figure consists of three panels labeled a, b, and c. All panels share a common Y-axis representing F I T in percent, ranging from 0 to 100 in increments of 10, and a common X-axis representing the Number of sensors, n sub p, ranging from 1 to 16. Three models are plotted in each panel: n 4 s i d plus P O D in red with square markers, linear 1 plus P O D in black with square markers, and linear 2 plus P O D in blue with square markers.

* Panel a, n equals r equals 20: The black and blue lines show a stable trend between 60 and 70 percent. The red line starts low at approximately 30 percent and shows a steep linear increase before plateauing around 65 percent at n sub p equals 11.

* Panel b, n equals r equals 25: The black and blue lines remain stable between 60 and 70 percent. The red line shows significant instability for low sensor counts, dropping to a trough of 10 percent at n sub p equals 4 before rising to stabilize near 55 percent.

* Panel c, n equals r equals 30: The black and blue lines continue to track closely between 60 and 70 percent. The red line exhibits fluctuating behavior between 30 and 55 percent, showing less stability than in panel a.

Figure 13 compares the reconstructed flow statistics with the true statistics. All three models capture the general shape of the spatial distribution of the Reynolds stresses. Closer inspection reveals that all three models reconstruct a slightly smaller region of high $ \left\langle {v}^{\prime }{v}^{\prime}\right\rangle $ and fluctuating kinetic energy $ k $ far downstream. They also underestimate the peak values of these two variables, but again the linear 2 model gives the most promising results. Better results can be obtained if models are trained at different values of $ \mathit{\operatorname{Re}} $ and then applied to a new (unseen) value.

Figure 13.

Contour plots of flow statistics at the off-design condition of $ \mathit{\operatorname{Re}}=800 $ : (left column) $ \left\langle {u}^{\prime }{u}^{\prime}\right\rangle $ (middle column) $ \left\langle {v}^{\prime }{v}^{\prime}\right\rangle $ , (right column) kinetic energy of the fluctuating field $ k $ . ( $ a $ – $ c $ , top row) DNS data, reconstruction using Kalman filter based on ( $ d $ – $ f $ , second row) the n4sid model with $ n=30 $ ( $ g $ – $ i $ , third row) linear 1 model and ( $ j $ – $ l $ , bottom row) linear 2 model with $ r=30 $ . All results are obtained with $ {n}_p=15 $ , only the leading five sensors are shown. Square (■) markers indicate the location of the $ {c}^{\prime } $ sensors.

Figure 13. Long description

The figure consists of 12 panels labeled a through l. The x-axis for all panels ranges from negative 2 to 10, and the y-axis ranges from 0 to 4. A gray rectangular obstacle is positioned at the origin in every plot.

* Columns: The left column (a, d, g, j) shows the Reynolds stress component u prime u prime. The middle column (b, e, h, k) shows v prime v prime. The right column (c, f, i, l) shows the kinetic energy k. Each column has a color scale bar at the top, ranging from dark blue (0) to dark red (0.8 for u prime u prime, 0.76 for v prime v prime, and 0.46 for k).

* Rows:

- Top row (a-c): D N S data showing smooth, continuous flow structures emanating from the obstacle.

- Second row (d-f): Reconstruction using the n 4 s i d model with n equals 30.

- Third row (g-i): Reconstruction using the linear 1 model.

- Bottom row (j-l): Reconstruction using the linear 2 model with r equals 30.

In rows two through four, five black square sensors are marked and numbered 1 through 5. Sensors 1 and 2 are located near the obstacle at x over h approximately 2.5, while sensors 3, 4, and 5 are clustered further downstream at x over h approximately 9. The reconstructed plots show high fidelity to the D N S data, capturing the primary wake structures and intensity peaks downstream of the obstacle.