Spatiotemporal Proper Orthogonal Decomposition of turbulent channel flow
Abstract
Proper Orthogonal Decomposition is applied to the wall layer of a turbulent channel flow () in a novel way, where empirical eigenfunctions are defined in both space and time. Due to the statistical symmetries of the flow, the eigenfunctions are associated with individual wavenumbers and frequencies. Selfsimilarity of the dominant eigenfunctions, consistent with wallattached structures, is established. We show that the most energetic modes are characterized by time scales in the range 200300 viscous wall units. The new decomposition provides a natural measure of the convection velocity of structures, with a characteristic value of 12 in the wall layer. Finally, we show that the energy budget can be split into specific contributions for each mode and provides a new closure formulation for the quadratic terms.
1 Introduction
Proper Orthogonal Decomposition (POD) was first introduced in turbulence by Lumley (1967). Its derivation stemmed from the KarhunenLoève (KL) decomposition (Loève, 1977) which represents a squareintegrable centered stochastic process in the time domain as an infinite linear combination of orthogonal functions. If is defined over a finite range, then
(1) 
where is stochastic and are continuous functions. It is important to note that represents a stochastic variable and not a sample. Realizations of will be noted . The functions are the eigenfunctions of the covariance function , where the operator refers to expectation with respect to the measure of . In the KarhunenLoève derivation, the variable corresponds to time, but it could indicate any other variable  such as space.
Lumley (1967) adapted the decomposition to Fluid Mechanics: the samples were constituted by flow realizations, and the ergodicity assumption was used to replace the covariance function corresponding to an ensemble average with one obtained by the time average so that the KL transform was generally applied to space. He considered the spatial autocorrelation tensor , where now and represent the deterministic variable (space) and represents the ensemble average (which is simply the time average here). therefore represents the spatial autocorrelation tensor at zero time lag. As pointed out by George (2017), it is important to realize that the spatial autocorrelation tensor is different from its sampled estimation , where is the number of samples. In general, exact eigenfunctions of cannot be computed since the can only be approximated. However, for spatially homogeneous flows such as the channel flow in horizontal directions (Aubry et al., 1988), exact solutions are known a priori since POD modes are Fourier modes in the homogeneous directions.
A key insight given by George (2017) is that if the process is stationary in time, the KL or POD modes are Fourier modes in the deterministic variable . Following George (2017), it appears that the proper way to apply Proper Orthogonal Decomposition would be to apply it in both time and space, reverting to ensemble average (as in the original definition) to evaluate the covariance function. This will allow us to identify modes that can be associated with individual spatial wavenumbers and temporal frequencies. Such decompositions could be useful to identify key features and instability mechanisms underlying the flow dynamics and eventually attempt to control them. One particularly promising approach for this is the resolvent analysis (McKeon, 2017), which is based on the NavierStokes equations. In contrast, the present approach is databased, and can therefore provide a complementary viewpoint.
In the present paper, we apply fully spatiotemporal Proper Orthogonal Decomposition to the turbulent channel flow. Wall turbulence is characterized by a variety of spatiotemporal scales interacting in a highly complex fashion (Robinson, 1991; Dennis, 2015). It is well known that the flow is characterized by an alternation of high and lowspeed streaks aligned with the flow, with a typical spacing in the spanwise direction and a length of wall units (Kim et al. (1971), Stanislas et al. (2008), Jiménez (2013)), which are units based on the fluid viscosity and wallfriction velocity , and will be denoted with a + throughout the paper. These streaks are associated with longitudinal vortical structures pushing lowspeed fluid upwards and bringing highspeed fluid downwards, which results in a strong contribution to the Reynolds stress (Kim et al., 1971). This contribution is highly intermittent in space and time, with ‘bursts’ of turbulence production. Eulerian measures of the bursts yield time scales on the order of 100 wall units (Blackwelder & Haritodinis, 1983). It is not clear how this time scale is related to the characteristics of the coherent structures, which are associated with a selfregenerating cycle of about 300 wall units (Jiménez, 2013).
Classic spatial POD has been previously used to investigate the boundary layer (see e.g. Aubry et al. (1988), Moin & Moser (1989), Podvin & Lumley (1998)). Spatial eigenfunctions are determined from secondorder statistics. The amplitudes of eigenfunctions are characterized by a mixture of time scales and can only be computed by projecting the full instantaneous field onto the spatial eigenfunctions. In contrast, the new decomposition directly provides spatiotemporal patterns.
In this paper we focus on a relatively moderate Reynolds number , based on the fluid viscosity, channel halfheight and friction velocity. At these moderate Reynolds numbers our focus will be on the wall layer as large scales such as those observed by Smits et al. (2011) are not present in the flow. The rest of the paper is organized in the following manner. Section 2 gives details of the methodology for the new variant of POD. Section 3 presents the POD results. The contributions to the turbulent kinetic energy equation associated with each mode are examined in Section 4, followed by a conclusion in Section 5 which summarizes the key observations and results.
2 Description of the procedure
2.1 Numerical database
We use the finitevolume code SUNFLUIDH to simulate the incompressible flow in a channel. Further details on the configuration and numerical scheme can be found in Podvin & Fraigneau (2014). Periodic boundary conditions are enforced in the streamwise () and the spanwise direction (), and the channel lengths along these directions are and , where is the channel halfheight. The total domain is discretized using a grid of size , with a uniform grid for the horizontal directions and and a hyperbolic tangentdistributed mesh for the wallnormal direction . The velocity components in the streamwise, wallnormal and spanwise direction will be denoted respectively by or equivalently by . The simulation is conducted at , based on the friction velocity and channel halfheight , which corresponds to the simulation of Moser et al. (1999). The database consists of 16500 flow realizations, separated by . The time spanned by the database represents about viscous time units.
In most of the paper, we focus on the domain , but we also considered the full boundary layer where . Due to the strong decrease of the energy spectra with the streamwise wavenumber, and given the typical extent of longitudinal streaks in the wall layer (about 6001000 wall units (Jiménez, 2013)), the streamwise extent of the domain was limited to representing 900 wall units. The full spanwise extent of the domain (about 1800 wall units) was considered.
The methodology is summarized in Figure 1. Since POD is also applied in time, the autocorrelation tensor needs to be computed from several independent realizations of the same experiment. An ideal sample would consist of the Fourier transform in space and time of the velocity field corresponding to different databases obtained at the same Reynolds number. In practice, owing to the cost of the simulation, we split a single database into several contiguous chunks of length , each of which constitutes a sample. The samples are therefore not independent realizations, but we assume that the period is sufficiently larger than the characteristic time scales of the flow, or at least sufficiently large to allow separation of the dominant time scales. Yet cannot be too large in order to allow for a reasonable number of samples to be constituted. We note that as in Moin & Moser (1989) and Podvin (2001), the number of samples is doubled by considering spanwise reflections and quadrupled by considering both lower and upper parts of the channel (including spanwise reflections).
Each sample then consists of a block of velocity fields uniformly sampled at a rate . For most of the test cases, we have fixed the value of to be 3000 wall units. Results are reported for four datasets extracted from the single database, which represent different domain sizes and sampling rates. The characteristics of the different configurations are given in Table 1. The sampling rate and number of samples were varied in in the configurations D1D3. These datasets correspond to the wall layer, . The fourth one (D4) corresponds to the full boundarylayer height, .
D1 
0.14  100  60  29.5  2950  
D2  0.14  100  120  29.5  2950  
D3  0.14  500  132  5.9  2950  
D4  1  500  132  5.9  2950 
2.2 Full spatiotemporal POD
Owing to the homogeneity of the statistics in the spatial directions and and their stationarity, POD modes are Fourier modes in the horizontal directions as well as in time. Here, for each configuration, the Fourier transform of the velocity is computed in the temporal and the homogeneous spatial directions for each sample corresponding to a set of fields. The th component of the velocity field () in the physical space is then written as
(2) 
where denotes the th velocity component in the Fourier space corresponding to the streamwise wavenumber , spanwise wavenumber and frequency . Proper Orthogonal Decomposition is then applied in the wallnormal direction for each triad :
(3) 
where, for any uplet ,
(4) 
Here, represents the usual ensemble average over the space of all possible flow realizations, is the Kroneckerdelta function, denotes a specific POD mode for each triad, and * denotes complex conjugation. Note that is independent of as the three velocity components are stacked into a single vector. The complex, stochastic coefficients are uncorrelated and their variance is equal to .
The eigenfunctions and eigenvalues are obtained by solving the following eigenproblem, where the autocorrelation is estimated from taking an ensemble average:
(5) 
Here, represents the single velocity vector with the three components stacked in it. Upon discretization in , represents the single eigenvector containing the three components. By construction, the eigenvectors are orthogonal with respect to the Euclidean inner product, and can be interpreted as the energy content in each mode. Similarly, for each triad , the modes (i.e., different ’s) are sorted by energy. Since the lefthandside of Equation (5) is a HilbertSchmidt integral operator, its eigenvalues are real and nonnegative. The entire procedure is schematically represented in Figure 1. The eigenvalues associated to the triplets can then be gathered and sorted globally according to their magnitude. We will call mode number the global index associated with the sorted modes over spatial wavenumbers , frequency and quantum number , with lower mode number indicating a larger eigenvalue. In all that follows, we focus only on the most energetic modes.
3 Results from spatiotemporal POD
3.1 POD eigenvalues
Figure 1(a) shows the top 5000 eigenvalues for each configuration. We note that this represents a tiny fraction of the total number of eigenvalues defined in the spatiotemporal space, which is , where is the number of instantaneous fields contained in a sample, and are respectively the numbers of streamwise and spanwise wavenumbers, and is the number of grid points in the wallnormal direction. For D3 this corresponds to about degrees of freedom. The largest eigenvalue where and corresponds to that of the mean mode (discussed in the next section). In the rest of the paper we will focus on the eigenfunctions. We note that the eigenfunctions corresponding to capture about of the turbulent kinetic energy, which is in agreement with Moin & Moser (1989)’s results.
Figure 2a shows that the spectra in D1 and D2, which correspond to a different number of samples, are essentially indistinguishable from each other. This indicates that the number of samples of ensemble averaging appears to be sufficient for the convergence of the dominant eigenvalues. The spectra in D1 and D2 are similar to that in D3 at low values of , but have higher energy levels at large values of , which is likely to be due to aliasing effects as the fields are sampled at higher rates there than in D3. The spectrum has an asymptotic decay of . Unlike the other configurations, D4 corresponds to the full boundary layer . As expected, energy levels are higher but the shape of the spectrum and its asymptotic decay rate are the same as in the wall region , which suggests selfsimilarity (see next section). Figure 2b shows the fraction of turbulent kinetic energy (TKE) captured by the first modes for each (as we are considering fluctuations, the contribution from the mean mode was set to zero). The fraction of TKE captured by the most energetic mode numbers is relatively high, but increases only slightly for mode numbers higher than , which highlights the complexity of the flow. This is confirmed by Table 2, which shows that the first eigenvalues capture about of the TKE which is increased only to when modes are included, for the D3 configuration. The trend remains the same for the D4 configuration, where eigenvalues capture around of TKE which increases only to when eigenvalues are included.
Number of eigenvalues ()  % of TKE (D3)  % of TKE (D4) 

200  
3000  
5000 
Figure 3 shows the dependence of the integrated spectra with respect to the frequency and spanwise wavenumber for the different configurations. For the D3 configuration, the spatial wavenumbers and the frequency , expressed in wall units, are related to the mode indices as
(6) 
Figure 3a shows how the sum of the eigenvalues over streamwise wavenumber and frequencies () varies as a function of . Results show a good agreement with the rescaled standard energy spectrum obtained by Moser et al. (1999) at . As expected, the spanwise spectrum over the full layer has more energy in the lower wavenumbers and less energy at higher wavenumbers compared to that in the wall layer.
Figure 3b shows the sum of the dominant eigenvalues over all modes (i.e., ) as a function of frequency for the different configurations D1D4. As observed previously, aliasing effects can be observed for the fields sampled in time at a lower rate (D1D2), but the trends are very similar. A marked increase in the energy is observed in the frequency range of , which corresponds to time scales of 200300 viscous units. This value is in good agreement with the duration of the regeneration cycle identified by Hamilton et al. (1995) and Jimenez & Moin (1991) in minimal domains, as well as the investigations of largerscale domains reported in Jiménez (2013). A similar frequency peak is observed in the full boundary layer (D4), but with a shift towards slightly higher frequencies: the maximum is located at , while it is located at , in the region .
The predominance of the characteristic frequency at is confirmed by Table 3, which shows the top 30 eigenvalues (denoted by ) with the corresponding Fourier mode indices in space and time and quantum mode for the configuration D3. Although a large fraction (18 modes) of the first 30 modes are associated to , which corresponds to long time scales that are outside the scope of our analysis, all the other modes are characterized by frequencies between and , which correspond to the previously identified time scale of about 200300 viscous units. The most energetic of these modes corresponds to , i.e., . We note that all the modes in the table are characterized by , which corresponds to the largest eigenvalue for a triad , and , which corresponds to the streamwiseaveraged flow.
1  0  1  24.97645  
2  0  1  0.008133  
3  0  1  0.004205  
4  0  1  0.003754  
5  0  1  0.001489  
6  0  1  0.001277  
7  0  1  0.001222  
8  0  1  0.001128  
9  0  1  0.000923  
10  0  1  0.000655  
11  0  1  0.000591  
12  0  1  0.000552  
13  0  1  0.000523  
14  0  1  0.000505  
15  0  1  0.000457 
16  0  1  0.0004568  
17  0  1  0.0004559  
18  0  1  0.000451  
19  0  1  0.000415  
20  0  1  0.000412  
21  0  1  0.000396  
22  0  1  0.0003806  
23  0  1  0.0003643  
24  0  1  0.000349  
25  0  1  0.0003475  
26  0  1  0.0003445  
27  0  1  0.000335  
28  0  1  0.000314  
29  0  1  0.000295  
30  0  1  0.000289 
To provide a comparison with standard analysis tools, Figures 4a and b show the temporal Fourier transform of the streamwise velocity component averaged along horizontal directions for two different heights: and . Both spectra show, although not distinctly, a peak around the characteristic frequency . This illustrates the usefulness of the new POD implementation.
Figures 5a and 5b respectively show the distribution of the top 5000 eigenvalues along and for the D3 configuration. Each eigenvalue is represented by a dot. Both figures show that the peaks observed in frequency and wavenumber space are not created by an accumulation of small eigenvalues, but correspond to coherent, more energetic structures. We note that Figure 5b also shows smaller, but noticeable peaks of at around 0.009 and 0.0135, which correspond to the harmonics of the frequency .
Figure 5(a) shows that the energy spectrum integrated in frequency space slowly decreases in an apparently selfsimilar manner with respect to both streamwise and spanwise wavenumbers. Figure 5(b) shows the corresponding variations in the frequency space of the energy spectrum integrated in the spanwise wavenumber space for different streamwise wavenumbers . For the mode, owing to spanwise reflection invariance, the plot is symmetric in the frequency space, and the characteristic frequency and its harmonics can be clearly identified. These peaks observed at are still present in the spectra, which confirms that the characteristic frequency is not an artefact of the streamwise average. However nonzero streamwise wavenumbers are also characterized by a broad peak, which represents convection effects. In general, defining a convection velocity for the turbulent wall layer is not straightforward (Krogstad et al. (1998); Alamo & Jimenez (2009)), but the question can be more easily addressed in the present framework where each POD eigenfunction is naturally associated with a phase velocity . The spectra can be used to define a global convection velocity at each streamwise wavenumber using
(7) 
where
(8) 
Peaks in the spectra are located at for respective streamwise wavenumbers of . This corresponds to a global convection velocity of 12, with a slight upward shift observed with increasing streamwise wavenumbers. This agrees well with the classical results (Kreplin & Eckelmann (1979), Krogstad et al. (1998), Wallace (2014), Alamo & Jimenez (2009)).
3.2 POD eigenfunctions
Figure 7 shows that the reconstructed mean streamwise profile coincides with previous results (Moser et al. (1999)) and with the dominant mode of the decomposition. Figure 8 shows the absolute value of each component of the most energetic mode , for , obtained for the datasets D1D3. As expected, there is a good agreement between the different configurations, which indicates convergence of the procedure at least for the most energetic modes. In all that follows only results for D3 and D4 will be presented.
In order to gain more insight on the structure of the modes, Figures 9a, c and e show the shape of the most energetic eigenfunction components associated with the characteristic frequency for different spanwise wavenumbers. The location of the velocity maxima moves closer to the wall as the spanwise wavenumber increases, which is in agreement with previous descriptions of wallattached structures (Alamo et al. (2006), Podvin et al. (2010)). The values of the maxima are similar for both crossstream components, also in agreement with previous observations (Alamo et al. (2006), Podvin & Fraigneau (2017)), but the value tends to increase with the spanwise wavenumber for the streamwise and the wallnormal component, while it slightly decreases for the spanwise component. The monotonous evolution of the shape of the modes with the spanwise wavenumber suggests selfsimilarity, as was proposed in Podvin & Fraigneau (2017).
Figures 9b, d and f compare the eigenfunctions obtained on the domain (for D3) and (for D4) for several spanwise wavenumbers. The eigenfunctions of D4 have been rescaled to have the same energy content as those of D3 on . We observe that the eigenfunctions nearly coincide over their common definition domain, which is not a trivial result. The persistence of the eigenfunction shape with respect to the wallnormal extension of the decomposition domain shows the coherence of the most energetic motions over the entire height of the boundary layer. It also means that the restrictions of the eigenfunctions on the larger domain are orthogonal to each other on the smaller domain (since they coincide with the eigenfunctions there), which makes it possible to recover the amplitude of the eigenfunction on the larger domain directly from the projection of the velocity field in the smaller domain onto the corresponding eigenfunction. This shows the relevance of the decomposition for estimation purposes in a context of partial information (see for instance Podvin et al. (2010)).
Figure 10 shows the flow generated by selected modes in figure 9 associated with for two different values of . As these dominant modes have no streamwise variation (i.e., ), the threedimensional shapes of the modes are represented by a contour plot of the streamwise velocity component () and a velocity plot of the inplane quantities . The modes are characterized by updrafts of lowspeed fluid alternating with downdrafts of highspeed fluid, associated with vortical motions in the streamwise direction, in agreement with classic observations (Corino & Brodkey, 1969). The wallnormal extension of the structure decreases as the spanwise wavenumber increases, in agreement with Townsend’s model (Townsend, 1976) of wallattached eddies (note the similarity of Figure 10 with Figure 9 in Jiménez (2013)). Figure 11 shows that the most energetic streamwise mode convected at a velocity of also corresponds to high and lowspeed streaks alternating in the streamwise as well as in the spanwise direction and associated with vortical motions in the crossstream plane. The wallnormal position of the vortex centers depends on the vertical extension of the streaks.
a) 
b) 
The contribution of the most energetic modes to the fluctuating stresses can be evaluated using
(9) 
where the subscript ’’ is used to distinguish it as a reconstructed quantity using most energetic modes. The bar indicates average over space and time and denotes the usual ensemble operator. Note that the mean modes corresponding to have been excluded from the reconstruction as the stresses are associated with fluctuating components. Figures 12a, b and c show the three components of the rootmean square velocity (rms) and Figure 12d shows the Reynolds stress as a function of reconstructed using 200 and then 3000 modes. The plots are compared with the present simulation results and the DNS results of Moser et al. (1999). A nonnegligible fraction of the turbulent intensity is recovered with the most energetic 200 modes: 60% for the streamwise velocity, 30% for the crossstream components, and about 25% of the Reynolds stresses. The gain obtained using 15 times as many modes (3000) is relatively small (40% increase for the spanwise component and 30% for the wallnormal component), which highlights the complexity of the flow.
4 Energetics of modes
With the new variant of Proper Orthogonal Decomposition, projection of the NavierStokes equations onto the basis of eigenfunctions no longer yields a dynamical system for the temporal amplitudes of the spatial eigenfunctions, but makes it possible to evaluate the contributions made by the different terms of the momentum equations for each mode .
For consistency with standard analysis, we decompose the th component of the velocity field into its mean part which, as we have seen above, can be represented by the mode and a fluctuating part .
We rewrite the equations as
(10) 
where (mean streamwise velocity) and if and only if and .
To obtain a budget for each mode we project Equation (10) onto the corresponding mode, i.e., we take the Fourier transform in the space, take the inner product in the wallnormal direction, and apply an ensemble average. This gives
(11) 
where denotes complex conjugation.

The first term yields a contribution . We note that this contribution is purely imaginary.

The second and third terms represent the interaction of the mode with the mean velocity profile. The second term is associated with the classic production term and is expected to represent a source of energy for the fluctuations:
(12) The third term represents the convection effect of the mean field:
(13) 
The last term on the lefthandside of the equation corresponds to the viscous diffusion term defined as:
(14) Integrating the last term by parts, one has
(15) Except for the last term corresponding to a boundary effect, all contributions to are real and negative, which corresponds to an energy loss for the mode, as expected.
Note that and only require information about the mode (one can think of them as purely diagonal operators) and can be directly evaluated from the Proper Orthogonal Decomposition, while the other two contributions require additional information about the coefficients .

As in the classic derivation (Aubry et al. (1988)) , the contribution from the pressure term represents the influence of the pressure at the upper boundary of the wall layer and can be expressed as
where is the Fourier transform (along , and ) of pressure at height , which needs to be evaluated from the DNS. This term is proportional to the wallnormal intensity of the structure at the top of the layer and depends on the velocitypressure correlation at the top of the layer. It represents an external forcing term which corresponds to the interaction of the wall layer with the outer region (Aubry et al., 1988).

Finally, the quadratic interactions can be evaluated as
(16) requires information about the triple correlations of the coefficients , and involves first derivatives of the eigenfunctions. This term characterizes the energy transfer from the different modes to the mode . It can be seen as a nonisotropic extension of the energy transfer function defined in isotropic turbulence (Zhou, 1993), and involves triads of modes in the space. It also corresponds to the forcing term of the resolvent analysis (McKeon, 2017).
If all other terms are known, the quadratic interaction terms can be determined from the budget equation:
(17) The real part of the equation represents a balance between the production and dissipation term , which depend only on the characteristics of the mode , and the interaction terms due to pressure and convection by velocity fluctuations , which characterize how the mode interacts with the full flow. The imaginary part of the equation can be seen as a phase dispersion relation linking the frequency of the mode () with the different physical mechanisms.
Figure 13 represents the contributions of the different terms to the equations for the largest modes (the contributions of the modes corresponding to the mean flow were set to zero). These modes are all characterized by a zero streamwise wavenumber () and a quantum number . All terms are evaluated directly, except the quadratic term which is evaluated using Equation (17). We have checked (not shown here though) that the quadratic terms could not be evaluated correctly by direct computation limited to the first 200 modes, as higherorder contributions were significant, which is typical of the closure problem.
For the sake of clarity, the left plots show mode number from and the right plots show from . The top row (Figures 13a and b) shows that the real part of the production term is essentially balanced by the sum of the real part of the dissipation and the quadratic terms. For the less energetic modes (Figure 13b) we note that the quadratic terms are mostly negative, which is consistent with the idea of a positive energy transfer from the large scales (most energetic modes) to the small scales. Figure 13 shows that the pressure interaction term is very small compared to the other terms, in agreement with Aubry et al. (1988)’s derivation where it is modelled as a stochastic term of small amplitude. If we neglect the influence of the pressure term we have
(18) 
Figures 13c and d show that for the most energetic modes, as , the frequency of the mode is directly related to the imaginary part (denoted by ) of the nonlinear contributions to the mode :
(19) 
This corresponds to the following closure approximation:
(20) 
As observed above, we can see that triple correlations between the different POD modes are essential to characterize the energy transfer between the scales. Computing these correlations is difficult, however if can be determined from Equation (17), the decomposition offers a new way to reconstruct the quadratic terms using
(21) 
Substitution of the velocity decomposition into the lefthand side of the equation requires performing a cumbersome convolution on all coefficients , but the righthandside provides a straightforward expression of the quadratic terms in the stochastic POD basis of coefficients , where coordinates are solely determined from secondorder statistics. At each scale level , the effect of the quadratic terms is to stretch and to rotate the corresponding velocity mode by a factor . Direct modelling of the distribution of could therefore lead to new formulations of turbulence models in the decomposition framework.
5 Conclusion
Spatiotemporal Proper Orthogonal Decomposition has been applied to the wall layer of a turbulent channel flow. The decomposition represents an efficient data reduction technique which is adapted to large simulation databases. It brings to light typical features of wall turbulence in a straighforward manner, but also provides a fresh viewpoint on the flow organization. Due to symmetry properties, the decomposition singles out empirical eigenfunctions for each frequency and horizontal spatial wavenumber. Besides time scales superior to 3000 wall units, which our limited implementation of POD did not allow us to characterize fully, we have shown that the most energetic modes were characterized by a time scale on the order of 250 wall units. About 30% of the turbulent kinetic energy was captured by the 200 most energetic modes. The most energetic modes were found to have a selfsimilar shape that appeared largely independent from the wallnormal extent of the decomposition domain, which shows the coherence of the motions over the height of the boundary layer. Convection velocities could be directly defined from the POD spectrum. A global convection velocity on the order of 12 was identified in the wall layer, in good agreement with previous approaches.
Finally, substitution of the decomposition into the NavierStokes equation and numerical computation of the different contributions of the modes to the turbulent kinetic energy budget highlighted the key role played by quadratic interactions and allowed us to propose a new closure formulation to model the contribution of these interactions. Such relationships, which need to be further explored in a careful manner, could be useful to derive new turbulence models. We hope that this work will pave the way for comprehensive investigations of wall turbulent flows using spatiotemporal Proper Orthogonal Decomposition.
Acknowledgements
This work was supported by the Center of Data Science from the ParisSaclay University. Computations were carried out at IDRISGENCI (project 02262).
References
 Alamo, J. C. Del & Jimenez, J. 2009 Estimation of turbulent convection velocities and corrections to taylor’s approximation. Journal of Fluid Mechanics 640.
 Alamo, J. C. Del, Jimenez, J., Zandonade, P. & Moser, R. D. 2006 Selfsimilar vortex clusters in the turbulent logarithmic region. Journal of Fluid Mechanics 561.
 Aubry, N., Holmes, P., Lumley, J. L. & Stone, E. 1988 The dynamics of coherent structures in the wall region of the wall boundary layer. Journal of Fluid Mechanics 192.
 Blackwelder, R. F. & Haritodinis, J. H. 1983 Scaling of the bursting frequency in turbulent boundary layers. Journal of Fluid Mechanics 132.
 Corino, E. R. & Brodkey, R. S. 1969 A visual investigation of the wall region in turbulent flow. Journal of Fluid Mechanics 37.
 Dennis, David J.C. 2015 Coherent structures in wallbounded turbulence. Anais da Academia Brasileira de Ciencias 87, 1161 – 1193.
 George, W. K. 2017 A 50year retrospective and the future. In Whither Turbulence and Big Data in the 21st Century?. Springer.
 Hamilton, J. M., Kim, J. & Waleffe, F. 1995 Regeneration mechanisms of nearwall turbulence structures. Journal of Fluid Mechanics 287.
 Jiménez, J. 2013 Nearwall turbulence. Physics of Fluids 25 (10).
 Jimenez, J. & Moin, P. 1991 The minimal flow unit in nearwall turbulence. Journal of Fluid Mechanics 225.
 Kim, H. T., Kline, S. J. & Reynolds, W. C. 1971 The production of turbulence near a smooth wall in a turbulent boundary layer. Journal of Fluid Mechanics 50 (1).
 Kreplin, H. P. & Eckelmann, H. 1979 Propagation of perturbations in the viscous sublayer and adjacent wall region. Journal of Fluid Mechanics 95.
 Krogstad, P. A., Kaspersen, J. H. & Rinestead, S. 1998 Convection velocities in a turbulent boundary layer. Physics of Fluids 10.
 Loève, M. 1977 Probability Theory. Berlin, Germany: Springer.
 Lumley, J. L. 1967 The structure of inhomogeneous turbulent flows. In Atmospheric Turbulence and Radio Wave Propagation. Nauka, Moscow.
 McKeon, B. J. 2017 The engine behind (wall) turbulence: perspectives on scale interactions. Journal of Fluid Mechanics 817.
 Moin, P. & Moser, R. D. 1989 Characteristiceddy decomposition of turbulence in a channel. Journal of Fluid Mechanics 200.
 Moser, R. D., Kim, J. & Mansour, N. N. 1999 Direct numerical simulation of turbulent channel flow up to = 590. Physics of Fluids 11 (4).
 Podvin, B. 2001 On the adequacy of the tendimensional model for the wall layer. Physics of Fluids 13 (1).
 Podvin, B. & Fraigneau, Y. 2014 PODbased wall boundary conditions for the numerical simulation of turbulent channel flows. Journal of Turbulence 15 (3).
 Podvin, B. & Fraigneau, Y. 2017 A few thoughts on proper orthogonal decomposition in turbulence. Physics of Fluids 29.
 Podvin, B., Fraigneau, Y., Jouanguy, J. & Laval, J. P. 2010 On selfsimilarity in the inner wall layer of a turbulent channel flow. Journal of Fluids Engineering 132 (4).
 Podvin, B. & Lumley, J. L. 1998 A lowdimensional approach for the minimal flow unit. Journal of Fluid Mechanics 362.
 Robinson, S. K. 1991 Coherent motions in the turbulent boundary layer. Annual Review of Fluid Mechanics 23.
 Smits, A. J., McKeon, B. J. & Marusic, I. 2011 High–reynolds number wall turbulence. Annual Review of Fluid Mechanics 43 (1).
 Stanislas, M., Perret, L. & Foucaut, J. M. 2008 Vortical structures in the turbulent boundary layer: a possible route to a universal representation. Journal of Fluid Mechanics 602.
 Townsend, A. A. 1976 The Structure of turbulent shear flow. Cambridge University Press.
 Wallace, J. M. 2014 Spacetime correlations in turbulent flow: A review. Theoretical and Applied Mechanics Letters 4.
 Zhou, J. 1993 Interacting scales and energy transfer in isotropic turbulence. Tech. Rep. CR191477. NASA.