1 Introduction
The ready availability of simulation and experimental data has led to a surge in interest in data-driven model reduction within the fluid mechanics community. Among their many potential uses, one of the common applications of these methods is to obtain a low-order reconstruction of the flow data. Methods such as proper orthogonal decomposition (POD) (Sirovich Reference Sirovich1987) and dynamic mode decomposition (DMD) (Schmid Reference Schmid2010) seek to accomplish this by identifying modes representing the dominant organized motions within the flow. The most broadly used form of POD generates spatially coherent modes that optimally capture the average flow energy measured in a spatial norm, while DMD generates temporally coherent modes, each associated with a single complex frequency, that correspond to the eigenvectors of the best-fit linear operator that advances the flow state from one time instant to the next.
Recently, Towne, Schmidt & Colonius (Reference Towne, Schmidt and Colonius2018) showed that a variant of POD (Lumley Reference Lumley1970), often called spectral proper orthogonal decomposition (SPOD), combines many of the advantages of both POD and DMD for statistically stationary flows: it produces modes that evolve coherently in space and time (i.e. literal coherent structures) and optimally captures the two-point space–time correlations and the flow energy in a space–time norm.
Despite their optimality, the most energetic SPOD modes typically capture a modest fraction of the total flow energy in broadband turbulent flows, and many modes may be required to accurately reconstruct the flow or its statistics. For example, Towne et al. (Reference Towne, Colonius, Jordan, Cavalieri and Brès2015) and Schmidt et al. (Reference Schmidt, Towne, Rigas, Colonius and Brès2018) found that the leading SPOD mode at each frequency captures only approximately 10 %–20 % of the flow energy in a turbulent jet, and that more than fifty modes are required at each frequency to capture 95 % of the energy. Similar observations have been made for other turbulent flows, including wakes, especially at high frequencies (Sanjose et al. Reference Sanjose, Towne, Jaiswal, Moreau, Lele and Mann2019; Symon, Sipp & McKeon Reference Symon, Sipp and McKeon2019)
In this paper, we propose to use Gabor modes to enrich the low-order SPOD flow reconstruction in order to recover truncated portions of the flow. Gabor modes are compact support wavepackets, with each mode carrying its own real valued wavevector, 
$\boldsymbol{k}$, a real valued spatial location, 
$\boldsymbol{x}$, and complex valued velocity vector, 
$\hat{\boldsymbol{u}}$. As shown by Ghate & Lele (Reference Ghate and Lele2017) and Ghate (Reference Ghate2018), representation of subfilter-scale turbulence using Gabor modes has three distinct advantages: (a) the modes are equivalent to a POD basis for quasi-homogeneous or locally homogeneous turbulence, and such a discrete modal representation can be obtained using a very small number of unique modes, (b) temporal evolution of the modes can be derived as a WKB-asymptotic approximation to subfilter-scale Navier–Stokes equations, and is equivalent to solving an ordinary differential equation (ODE) for each mode travelling in a Lagrangian frame of motion, and (c) the transform from modal/wavespace to physical space can be performed efficiently using non-uniform fast Fourier transforms (NUFFTs). Critically, the Gabor modes used to enrich the SPOD reconstruction are a function only of the retained SPOD modes, and therefore do not increase the order of the reconstruction.
In summary, the method proposed herein enables flow field reconstruction using fewer SPOD modes. We emphasize that the approach can be applied to reconstruct incompressible, inhomogeneous, temporally stationary turbulence in arbitrary three-dimensional domains; however, the present application addresses the problem of reconstruction on a two-dimensional plane, due to its wider appeal in applications that require inflow generation, ranging from scale-resolving simulations in aerospace engineering (Deck Reference Deck2005; Brès et al. Reference Brès, Jordan, Jaunet, Le Rallic, Cavalieri, Towne, Lele, Colonius and Schmidt2018), to meso-microscale coupling problems in geophysics and wind energy (Muñoz-Esparza et al. Reference Muñoz-Esparza, Kosović, Mirocha and van Beeck2014).
2 Problem set-up
The problem being considered is that of interaction between broadband isotropic turbulence and a dragging actuator disk. This set-up was originally studied by Ghate et al. (Reference Ghate, Ghaisas, Lele and Towne2018) and is shown in figure 1. The fixed thrust coefficient dragging actuator disk corresponds to the model of Calaf, Meneveau & Meyers (Reference Calaf, Meneveau and Meyers2010), and the Case 10 HIT inflow from Ghate et al. (Reference Ghate, Ghaisas, Lele and Towne2018), which corresponds to an integral length scale of roughly 
$25\,\%$ of the actuator-disk diameter, is the focus of this paper. The problem is formulated in a periodic domain (enabling the use of a 
$2/3$rd de-aliased Fourier discretization) with the inflow generated via a forcing fringe region (Nordström, Nordin & Henningson Reference Nordström, Nordin and Henningson1999) that uses a concurrent simulation of forced homogeneous isotropic turbulence. The simulation is performed in the infinite Reynolds number limit (
$\unicode[STIX]{x1D708}\rightarrow 0$) with subgrid closure provided by the Sigma model (Nicoud et al. Reference Nicoud, Toda, Cabrit, Bose and Lee2011). The analysis and modelling presented in the remainder of this paper focuses on data extracted on a single transverse ((
$y$–
$z$)-plane) located 
$6.66D$ units downstream of the actuator disk.
Problem configuration. Inflow homogeneous isotropic turbulence is generated using a concurrent forced HIT simulation with the desired integral length scale and dissipation rate.
3 Modal representation of wake turbulence
The present work addresses the reconstruction problem for spatially inhomogeneous but temporally stationary turbulence by proposing a reconstruction of the turbulent flow field on an arbitrary (
$y{-}z$)-plane as
$$\begin{eqnarray}\boldsymbol{u}(\,y,z,t)=\overline{\boldsymbol{U}}(\,y,z)+\boldsymbol{u}^{SPOD}(\,y,z,t)+\boldsymbol{u}^{Gab}(\,y,z,t),\end{eqnarray}$$ where 
$\overline{\boldsymbol{U}}$ is the mean velocity (not the focus of present work), while 
$\boldsymbol{u}^{SPOD}$ and 
$\boldsymbol{u}^{Gab}$ are portions of the velocity field that are represented using SPOD and Gabor modes, respectively.
3.1 Spectral proper orthogonal decomposition
The flow field sampled on the (
$y{-}z$)-plane located 
$6.66D$ downstream from the actuator disk is first decomposed into an ordered set of orthogonal SPOD modes in polar coordinates (
$r{-}\unicode[STIX]{x1D703}$):
$$\begin{eqnarray}\hat{\boldsymbol{u}}(r,m,f)=\frac{1}{N_{T}M_{\unicode[STIX]{x1D703}}}\mathop{\sum }_{n=0}^{N_{T}}\mathop{\sum }_{j=0}^{M_{\unicode[STIX]{x1D703}}}\boldsymbol{u}(r,\unicode[STIX]{x1D703}_{j},t_{n})\text{e}^{\text{i}(m\unicode[STIX]{x1D703}_{j}+2\unicode[STIX]{x03C0}ft_{n}/T)}=\mathop{\sum }_{j=1}^{J}a_{j}(m,f)\unicode[STIX]{x1D733}_{j}(r,m,f),\end{eqnarray}$$ where 
$\unicode[STIX]{x1D733}_{j}(r,m,f)$ is the SPOD mode shape (in 
$r$) corresponding to the 
$m$th azimuthal wavenumber at a discrete frequency, 
$f$. The total number of modes, 
$J$, at each frequency–wavenumber pair is controlled by the number of realizations used to compute the SPOD modes, with each realization defined over a finite sampling interval, 
$T=5.12D/U_{\infty }$ discretized uniformly using 
$N_{T}$ sampling points (time steps). We take 
$J=117$ and 
$M_{\unicode[STIX]{x1D703}}=768$ (to avoid aliasing) in what follows; our numerical experiments with different choices of 
$J$ and 
$T$ showed that the most energetic wavepackets associated with the Kelvin–Helmholtz instability were adequately captured using this sampling. The modes are sorted according to their modal energy, 
$\unicode[STIX]{x1D706}_{j}(m,f)=\langle a_{j}^{\star }(m,f)a_{j}(m,f)\rangle$, which defines the contribution of each mode to the total kinetic energy of the flow. The modal energies and mode shapes are shown in figure 2 up to a Strouhal number of 
$25$. The Strouhal number corresponding to the Nyquist frequency of the temporal sampling is approximately 
$50$; the upper half of the frequency range is excluded for the present analysis to avoid spurious artefacts associated with numerical/discretization error and spatial de-aliasing.
Modal energies and shapes for SPOD of wake turbulence taken at a transverse plane located 
$6.66D$ downstream for the drag disk.
Figure 2(a), which shows the modal energies as a function of Strouhal number, suggests that a truncated representation consisting of the first 10 leading modes (and all values of 
$m$) is only likely to produce accurate second-order correlations for 
$St<1$. Furthermore, at the axial downstream location being considered (
$6.66D$), we do not observe any dominant tone (frequency) in the primary varicose mode (
$m=0$). Figure 2(b) suggests that there is substantially more energy in the 
$m=1$ mode compared to the 
$m=0$ mode, especially at low Strouhal numbers (
$St<1$). It is also interesting to see that high wavenumbers do contain a substantial amount of energy, as shown in figure 2(b); the energy decays only as a power law, as a function of the azimuthal wavenumber, 
$m$ (see figure 2b). The mode shapes shown in figure 2(c) suggest that the bulk of the energy in the low-order modes (small 
$m$ and 
$j$) at lower frequencies corresponds to the shear layer turbulence. The azimuthal homogeneity embedded by the Fourier representation in 
$\unicode[STIX]{x1D703}$ is very efficient at isolating the shear layer/inflectional turbulence from the free stream, including scales entrained and subsequently distorted by the mean shear in the wake. Finally, an important consequence of finiteness of the data available to compute this modal representation is the uncertainty associated with higher-order modes, since the overall representation only converges as 
$\sqrt{n_{samples}}$ (Welch Reference Welch1967), where 
$n_{samples}$ can only be increased by running longer simulations. As such, while higher-order (in both 
$m$ and 
$j$) modes do contain non-negligible amounts of energy, in order to estimate them, we need to run very long simulations on very-high-resolution numerical grids. In contrast, the lower-order modes at low Strouhal numbers can be obtained using few 
$n_{samples}$ using a sufficiently accurate coarse grid simulation with a good subgrid-scale closure.
3.2 Truncated SPOD
A truncated/filtered representation corresponds to the following expansion:
$$\begin{eqnarray}\boldsymbol{u}^{trunc}(\,y,z,t)=\boldsymbol{{\mathcal{I}}}_{(r,\unicode[STIX]{x1D703})}^{(\,y,z)}\left\{\mathop{\sum }_{|\,f|<f_{max}}\mathop{\sum }_{|m|<m_{max}}\mathop{\sum }_{j<j_{max}}a_{j}(m,f)\unicode[STIX]{x1D733}_{j}(m,f,r)\text{e}^{-\text{i}(m\unicode[STIX]{x1D703}+ft)}\!\right\},\end{eqnarray}$$ where 
$a_{j}(m,f)$ can be computed using the appropriate orthogonality relations (Towne et al. Reference Towne, Schmidt and Colonius2018), and 
$\boldsymbol{{\mathcal{I}}}_{(r,\unicode[STIX]{x1D703})\rightarrow (\,y,z)}$ is the interpolation operator. We will take velocity 
$\boldsymbol{u}^{SPOD}(\,y,z,t)$ in (3.1) to be the truncated SPOD expansion, 
$\boldsymbol{u}^{trunc}(\,y,z,t)$.
Decomposition of instantaneous flow (evaluated on a 
$180\times 180$ grid) into a truncated SPOD representation (evaluated on a 
$36\times 36$ grid) and the resulting residual fields (evaluated on a 
$180\times 180$ grid) at an arbitrary sampling time. The velocity components in the Cartesian frame are obtained using an SPOD expansion truncated at 
$m_{max}=30$, 
$j_{max}=15$ and 
$f_{max}=1.4$.
The truncated description being considered in the remainder of the paper uses 
$f_{max}=1.4$ (Strouhal number), 
$m_{max}=30$ (azimuthal wavenumbers) and 
$j_{max}=15$ (leading modes); however, we note that all arguments presented here apply for arbitrary choices of truncation parameters as long as the large-scale coherent motions, in this case related to the Kelvin–Helmholtz instability, are captured by the truncated SPOD expansion.
We can define a residual field, 
$\boldsymbol{u}^{res}$ as
$$\begin{eqnarray}\boldsymbol{u}^{res}(\,y,z,t)=\boldsymbol{u}(\,y,z,t)-\boldsymbol{u}^{trunc}(\,y,z,t),\end{eqnarray}$$ where 
$\boldsymbol{u}$ is the instantaneous (fluctuating) component in the independent sample. Due to the orthogonality properties of the SPOD representation, it is easy to show that the total domain-averaged Reynolds stresses are given as the following superposition:
$$\begin{eqnarray}\langle u_{i}u_{j}\rangle =\langle u_{i}^{trunc}u_{j}^{trunc}\rangle +\langle u_{i}^{res}u_{j}^{res}\rangle .\end{eqnarray}$$ Figure 3 shows an example of such a representation using an arbitrary sample from the simulation that was not used to compute the modes, 
$\unicode[STIX]{x1D733}_{j}(m,f)$. It is important to note that the truncated representation can be evaluated on a 
$32\times 36\times 36$ Cartesian grid in 
$(\,y{-}z{-}t)$ space without any aliasing, as opposed to a 
$180\times 180\times 512$ grid needed for the full fluctuating field. These results indicate that the expansion given in (3.3) serves as an excellent surrogate to isolate large-scale, space–time coherent flow features in the 
$y{-}z{-}t$ domain being considered. Figure 3 clearly suggests that the residual scales are primarily fine-scale features that also appear to display quasi-homogeneity – i.e. spatial homogeneity at length scales corresponding to the filtering length scale implied by the SPOD truncation. The profile of the residual single-point correlations (figure 6) further indicates that these small scales are devoid of major radial inhomogeneity, and hence can be interpreted as the scales corresponding to the distorted free stream turbulence.
4 Gabor mode enrichment
4.1 Stationary Gabor modes
In the present application, the temporal evolution equations (see (2.16)–(2.22) in Ghate & Lele (Reference Ghate and Lele2017)) for each Gabor mode can be simplified substantially by neglecting the inter-scale sweeping in the planar directions (
$y$ and 
$z$) since 
$(U_{y}+u_{y}^{SPOD})/(U_{x}+u_{x}^{SPOD}),(U_{z}+u_{z}^{SPOD})/(U_{x}+u_{x}^{SPOD})\ll 1$. The planar reconstruction region (
$y{-}z$) is decomposed into 
$18\times 18$ quasi-homogeneous regions, each seeded with 80 Gabor modes. Under these assumptions, each Gabor mode is simply assumed to be advected in the streamwise direction according to the local streamwise time-averaged velocity (Taylor’s hypothesis); hence we refer to these modes as stationary Gabor modes. The energy exchange between the mean and SPOD scales and Gabor modes is captured via the straining/distortion effect. The temporal evolution of each Gabor mode located at 
$(\,y,z)$ carrying a complex valued velocity, 
$\hat{\boldsymbol{u}}$, and a real valued wavevector, 
$\boldsymbol{k}$, from time step 
$N$ to 
$N+1$ separated by 
$\unicode[STIX]{x1D6E5}_{t}$ can be summarized by the following four-step procedure (see Ghate & Lele Reference Ghate and Lele2017 for further details):
$$\begin{eqnarray}\left.\begin{array}{@{}c@{}}\displaystyle \hat{u_{i}}^{\star }=\exp (\text{i}k_{x}\overline{U}_{x}(\,y,z)\unicode[STIX]{x1D6E5}_{t})\hat{u_{i}}^{N},\\ \displaystyle \hat{u} _{i}^{\star \star }=\left[\unicode[STIX]{x1D6FF}_{ij}+\unicode[STIX]{x1D6E5}_{t}\left(\left(\frac{2k_{i}^{N}k_{m}^{N}}{k_{p}^{N}k_{p}^{N}}-\unicode[STIX]{x1D6FF}_{im}\right)\frac{\unicode[STIX]{x2202}U_{m}}{\unicode[STIX]{x2202}x_{j}}(\,y,z,t)-\unicode[STIX]{x1D708}_{t}(k)(k_{p}^{N}k_{p}^{N})\unicode[STIX]{x1D6FF}_{ij}\right)\right]\hat{u_{j}}^{\star },\\ \displaystyle k_{i}^{N+1}=\left[\unicode[STIX]{x1D6FF}_{ij}-\unicode[STIX]{x1D6E5}_{t}\frac{\unicode[STIX]{x2202}U_{j}}{\unicode[STIX]{x2202}x_{i}}(\,y,z,t)\right]k_{j}^{N},\\ \displaystyle \hat{u} _{i}^{N+1}=\left[\frac{k_{i}^{\,N+1}k_{j}^{\,N+1}}{k_{m}^{N+1}k_{m}^{N+1}}-\unicode[STIX]{x1D6FF}_{ij}\right]\hat{u} _{j}^{\star \star },\end{array}\right\}\end{eqnarray}$$ where 
$U_{m}(\,y,z,t)=\overline{U}_{m}+u_{m}^{SPOD}(\,y,z,t)$. In (4.1), the first stage corresponds to advection of enriched turbulence in the direction normal to the sampling inflow plane (
$y{-}z$) by the time-averaged velocity. The second and third steps represent the straining of enriching small scales by the larger SPOD (and mean) scales and the modification of the wavevector is a forward Euler approximation to the Eikonal equation. Finally, the projection implied in the fourth step is primarily used to discretely impose the divergence-free constraint; since the second and third steps are forward Euler approximations to the governing ODEs (see Ghate & Lele Reference Ghate and Lele2017), they inherently possess a spurious divergence (
$O(\unicode[STIX]{x1D6E5}_{t}^{2})$), which can be removed at virtually no additional computational cost. The time step is chosen based on the smallest scale enriched (
$k_{max}$) and the advective velocity 
$\overline{U}_{x}$. The choice of 
$k_{max}$ is rather arbitrary, and based on the Nyquist criterion of the physical space numerical grid on which the enriched fields are rendered.
To initialize the Gabor modes in each quasi-homogeneous region, we begin by randomly sampling isotropic modes over log-spaced wavenumber shells with a prescribed energy spectrum (see (2.26a) in Ghate & Lele (Reference Ghate and Lele2017)) parameterized using a dissipation rate, 
$\unicode[STIX]{x1D700}$, and a length scale measure, 
$L_{iso}$; both parameters vary only radially in the present application. The dissipation rate is modelled as
$$\begin{eqnarray}\unicode[STIX]{x1D700}(r)=\unicode[STIX]{x1D700}_{\infty }-\langle u_{i}^{SPOD}u_{j}^{SPOD}\rangle \frac{\unicode[STIX]{x2202}\langle U_{i}\rangle }{\unicode[STIX]{x2202}x_{j}},\end{eqnarray}$$ where 
$\unicode[STIX]{x1D700}_{\infty }$ is the dissipation rate of the turbulent co-flow (ambient/free stream), which is typically known or can be computed using a RANS model or an SGS model. The length scale measure is computed as
$$\begin{eqnarray}L_{iso}(r)=c_{L}\unicode[STIX]{x1D70F}(r)[\langle U\rangle (r)],\end{eqnarray}$$ where the constant 
$c_{L}=1/0.816$ ensures that 
$L_{iso}(r\rightarrow \infty )$ corresponds to the integral length scale of the isotropic co-flow. The integral time scale, 
$\unicode[STIX]{x1D70F}$, is computed as the integral time scale of the large-scale axial velocity, 
$u_{x}^{SPOD}$. Figure 4 shows the profiles of the two model inputs as computed using the SPOD data. These isotropic modes are then distorted using the local mean velocity profile in accordance to rapid distortion theory through a wavenumber-dependent time scale (Mann Reference Mann1994; Ghate & Lele Reference Ghate and Lele2017), which results in anisotropic, small-scale turbulence that is consistent with the mean velocity gradients in the quasi-homogeneous regions. An example description of this procedure is shown in Ghate & Lele (Reference Ghate and Lele2017) in the context of sheared boundary layer turbulence.
Model inputs used to generate stationary Gabor modes. Dissipation is normalized as 
$\unicode[STIX]{x1D700}(r)/\unicode[STIX]{x1D700}(r\rightarrow \infty )$ and the length scale measure is normalized as 
$L_{iso}(r)/L_{iso}(r\rightarrow \infty )$.
Once initialized, the dynamics represented by (4.1) account for the following physical processes: (a) rapid time-scale energy transfer from the large SPOD scales into the enriched small scales that occurs due to large-scale strain, (b) consistent temporal decorrelation of small scales, since this occurs primarily due to large scales sweeping enriched scales, (c) effect of pressure as a Lagrange multiplier to impose the divergence-free constraint (ensured by the Eikonal equation for 
$\boldsymbol{k}$), and (d) decay of the intense small-scale Burgers vortices generated by the non-local (in scale space) interactions of the straining term by representing the local interactions (in scale space) due to the nonlinear relaxation as a spectral viscosity obtained using a renormalization group (RNG) model (Canuto & Dubovikov Reference Canuto and Dubovikov1996).
Finally, it is important to emphasize the computational efficiency of the enrichment algorithm. The overall computational cost can be decomposed into two steps: (a) temporal evolution of Gabor modes and (b) rendering (transform into physical space). Log-spaced sampling of Gabor modes results in substantial compression in representing small-scale turbulence (
${>}95\,\%$ in three dimensions and 
${>}20\,\%$ in two dimensions) and the time advancement for each mode is an entirely local operation as detailed in (4.1).
The rendering step which is required to obtain the enriched velocity field on a numerical mesh requires a NUFFT; in the present application the cost of each two-dimensional transform is equivalent to approximately 5–6 uniform two-dimensional FFTs. Further details of the algorithm are provided in Ghate (Reference Ghate2018).
While the application discussed in this paper focuses on planar reconstruction of wake turbulence, the use of Gabor modes for enrichment is more broadly applicable to a variety of three-dimensional complex flows, including wall-bounded turbulence. The algorithm requires some representation of domain/geometry-influenced large scales, either via a coarse large eddy simulation (LES) (Quon, Ghate & Lele Reference Quon, Ghate and Lele2018) or other data-driven techniques such as deep neural networks (Srinivasan et al. Reference Srinivasan, Guastoni, Azizpour, Schlatter and Vinuesa2019), which explicitly influence the small-scale dynamics modelled by the Gabor mode representation of the flow. The choice of SPOD basis used in the present work to represent temporally stationary large-scale flow physics is particularly convenient due to its orthogonality properties and spectrally sharp time-filtering of the truncated expansion.
4.2 Enrichment
Figure 5 shows an instantaneous snapshot of the inflow field taken at the same time as the one shown in figure 3. A qualitative comparison of the Gabor-mode-induced fields with the residual fields suggests good overall agreement. Perhaps the one striking differentiating feature is the somewhat higher azimuthal imprinting in the Gabor-mode-induced instantaneous fields compared to the instantaneous residual fields (see figure 3). While the small scales induced by the Gabor modes are indeed coupled with the truncated SPOD fields, due to the localized straining (stage 2 in (4.1)), the overall azimuthal symmetry is a consequence of the model inputs, which only vary radially.
Instantaneous snapshot of enriched flow fields. Sample time is comparable with that in the representative snapshots shown in figure 3.
Single-point correlations for the enriched fields.
Temporal auto-spectra for the three velocity components extracted at two radial locations of 
$r/R=1$ and 
$1.5$, where 
$R$ is the radius of the wake-generating actuator disk.
Azimuthal auto-spectra for the three velocity components extracted at two radial locations of 
$r/R=1$ and 
$1.5$, where 
$R$ is the radius of the wake-generating actuator disk.
In order to facilitate more quantitative comparisons of the enriched field with the true full fields, several statistical measures are shown in figure 6 through figure 8. Figure 6 shows the single-point correlations computed for the Gabor-mode-induced and residual fields, time- and ensemble-averaged. While the number of ensembles used for the residual fields is not sufficiently large to fully converge the statistics, the primary purpose of this figure is to demonstrate that an 
$18\times 18$ grid of Cartesian quasi-homogeneous regions on the (
$y{-}z$)-plane is sufficient to obtain the expected azimuthal symmetry in statistics. Fewer quasi-homogeneous regions would result in reduced spatial localization for the induced flow fields.
Assuming azimuthal symmetry, single-point correlations as a function of radial location are shown in figure 7. While the truncated SPOD expansion significantly underpredicts the correlations, the Gabor-mode-enriched field shows excellent statistical agreement with the original field. Through the bulk of the shear layer, the physical anisotropy (
$\langle u_{x}u_{r}\rangle$) is well captured by the Gabor modes; the slight underprediction of the turbulent kinetic energy at the shear layer centreline is notable. Upon closer inspection of the model inputs, we can explain this deficiency in terms of the estimated dissipation rate in (4.2). This definition predicts 
$\unicode[STIX]{x1D700}(r=0)=\unicode[STIX]{x1D700}(r\rightarrow \infty )$, which appears to be a substantial underprediction. In this model, the core turbulence entrained by the shear layer is not neglected. Since 
$E(k)\propto \unicode[STIX]{x1D700}^{2/3}$, any underprediction of the dissipation rate results in underprediction of variances.
The one-dimensional power spectra in frequency (Strouhal number) and azimuthal wavenumber (
$m$) are shown for the velocity components in figures 8 and 9 at two different radial locations. These spectra further corroborate the effectiveness of the current approach. It is interesting to note from figure 9 that the enrichment using Gabor modes leads to an increase in energy for 
$m<30$, since it is able to generate smaller azimuthal scales that are truncated due to filtering in time and the radial direction in the truncated expansion. Further discussion of one-dimensional spatial spectra in the context of subfilter-scale enrichment can be found in Ghate (Reference Ghate2018) for highly anisotropic near-wall turbulence.
Overall, our results show that spatially inhomogeneous turbulent flows can be effectively reconstructed by combining a few SPOD modes to capture the energy-containing coherent modes, which capture the large-scale inhomogeneity, enabling enrichment using Gabor modes. While, in Ghate & Lele (Reference Ghate and Lele2017), Gabor mode enrichment was assessed on wall-bounded turbulent flows using filtered LES data, it is promising to note the ability of the algorithm to accurately provide enrichment for a more conventional data-driven reduced-order modelling algorithm. A potential step towards future improvement is to address the energy deficiency that is seen in azimuthal and radial velocity components, near the cutoff frequency of the truncated SPOD reconstruction (
$St=1.4$). This is a consequence of an inconsistency between the geometric anisotropy implied by the resulting quasi-homogeneous regions (parameterized by 
$F_{co}$) and the true Reynolds-stress anisotropy of the subfilter scales. For the present choice of 
$F_{co}=1.4$, the resulting aspect ratio of the 
$[t,y,z]$ quasi-homogeneous regions is rather skewed (
${\approx}[5\times 1\times 1]$); this can be mitigated by a Reynolds-stress-informed choice of 
$F_{co}$.
5 Conclusions
A flow reconstruction method that combines data-driven modal analysis with physics-based turbulence enrichment is developed and tested for incompressible wake turbulence. For the actuator-disk wake considered in this paper, the circumferential symmetry of the shear layer is leveraged to represent the shear-layer-driven turbulence using a compressed set of SPOD modes. The orthogonality of SPOD modes allows us to interpret such a truncated representation as a filtering operation, which subsequently enables generation of subfilter scales via Gabor mode enrichment. This juxtaposition of data-driven modelling with physics-based enrichment enables efficient representation of statistically stationary flow fields that contain both large-scale coherent motions associated with inflectional instabilities and broadband 
$k^{-5/3}$ turbulence. Hence, the coupled formulation is more broadly applicable to a variety of statistically stationary turbulent flows including wall-bounded turbulence. We further emphasize that once the SPOD mode shapes are determined using data, an arbitrarily large number of random ensembles of statistically equivalent flow realizations can be generated via randomizing the phase of the complex valued amplitudes, 
$a_{j}(m,f)$, in (3.3). Each of these random realizations can be further enriched on-the-fly with smaller scales using randomly sampled Gabor modes, which provide a consistent extrapolation of the spectral content present in the SPOD representation. This procedure could be useful for generating ensembles of statistically equivalent inflow conditions containing both inhomogeneous large-scale and homogeneous small-scale motions.
Acknowledgements
A.S.G. was funded by Tomkat Center for Sustainable Energy at Stanford University. S.K.L. acknowledges partial support from NSF-CBET-1803378. All simulations were performed on the Stampede2 supercomputer under the XSEDE project ATM170028. The authors would also like to thank the anonymous referees for their thoughtful comments and contribution to this work.
Declaration of interests
The authors report no conflict of interest.
