Srikanth Derebail Muralidhar

Bérengère Podvin

Lionel Mathelin

Yann Fraigneau

Spatio-temporal Proper Orthogonal Decomposition of turbulent channel flow

An extension of Proper Orthogonal Decomposition is applied to the wall layer of a turbulent channel flow (Re τ = 590), so that 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. Self-similarity of the dominant eigenfunctions, consistent with wall-attached structures transferring energy into the core region, is established. The most energetic modes are characterized by a fundamental time scale in the range 200-300 viscous wall units. The full spatio-temporal decomposition provides a natural measure of the convection velocity of structures, with a characteristic value of 12u τ in the wall layer. Finally, we show that the energy budget can be split into specific contributions for each mode, which provides a closed-form expression for nonlinear effects.

Introduction

Proper Orthogonal Decomposition (POD) was first introduced in turbulence by [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF]. Its derivation stemmed from the Karhunen-Loève (KL) decomposition [START_REF] Loève | Probability Theory[END_REF] which represents a square-integrable centered stochastic process in the time domain U (t) as an infinite linear combination of orthogonal functions. If t is defined over a finite range, then

U (t) = n a n χ n (t), (1.1)
where a n is stochastic and χ n are orthogonal, square integrable, functions. In all that follows the superscript refers to the mode index. It is important to note that U (t) represents a stochastic variable and not a sample. Realizations of U (t) will be noted u(t).

The functions χ n (t) are the eigenfunctions of the covariance function

K U (t, t ) = E[U t U t ],
where the operator E refers to expectation with respect to the measure of U . In the Karhunen-Loève derivation, the variable t corresponds to time, but it could indicate any other variable -such as space. [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF] (see also [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF]) 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 † Email address for correspondence: Berengere.Podvin@limsi.fr arXiv:1805.01494v3 [physics.flu-dyn] 9 Jan 2019 the one obtained by the time average so that the KL transform was generally applied to space. He considered the spatial autocorrelation tensor K U (x, x ) = U (x, t)U (x , t) , where now x and x represent the deterministic variable (space) and • represents the ensemble average (which is simply the time average here). K U (x, x ) therefore represents the spatial autocorrelation tensor at zero time lag. As pointed out by George (2017), it is important to realize that the spatial auto-correlation tensor K U (x, x ) is different from its sampled estimation 1 Ns Ns i=1 u(x, t i )u(x , t i ), where N s is the number of samples. In general, exact eigenfunctions of K U (x, x ) cannot be computed since the K U (x, x ) can only be approximated. However, for spatially homogeneous flows such as the channel flow in horizontal directions [START_REF] Aubry | The dynamics of coherent structures in the wall region of the wall boundary layer[END_REF], exact solutions are known a priori since POD modes are Fourier modes in the homogeneous directions.

A key insight stated by George ( 2017) is that if the process is stationary in time, the KL or POD modes are Fourier modes in the deterministic variable t. Combining Fourier transform in time with Proper Orthogonal Decomposition in space was first applied in experimental studies of free turbulent shear flows, such as mixing layers and jets. Spatial and frequency decomposition was performed in the pioneering work of [START_REF] Glauser | An application of lumley's orthogonal decomposition to the axisymmetric jet mixing layer[END_REF], [START_REF] Glauser | An orthogonal decomposition of the axisymmetric jet mixing layer utilizing cross-wire measurements[END_REF] as well as in [START_REF] Arndt | The proper orthogonal decomposition of pressure surrounding a turbulent jet[END_REF], [START_REF] Delville | Examination of large-scale structures in a turbulent plane mixing layer. part 1. proper orthogonal decomposition[END_REF], [START_REF] Citriniti | Reconstruction of the global velocity field in the axisymmetric mixing layer utilizing the proper orthogonal decomposition[END_REF], [START_REF] Ukeiley | Examination of large-scale structures in a turbulent plane mixing layer. part 2. dynamical systems model[END_REF] to cite only a few. However, as far as we know, a full four-dimensional decomposition was never attempted for any flow until recently (see the spectral POD of [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]), and has never been implemented for wall-bounded flows. The present work builds on these previous developments, and implements George's suggestion to apply the decomposition in the four-dimensional space, i.e. 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 [START_REF] Mckeon | The engine behind (wall) turbulence: perspectives on scale interactions[END_REF], which is based on the Navier-Stokes equations. In contrast, the present approach is data-based, and can therefore provide a complementary viewpoint.

As mentioned above, we note that another four-dimensional decomposition is provided by the Spectral POD introduced by [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. Their method only requires stationarity in time, so that temporal Fourier transform can be applied, while standard snapshot POD is performed in the spatial directions. In contrast, in our implementation, we apply Fourier transform in the spatial (horizontal) directions where turbulence is homogeneous. The connection of the present method with the snapshot POD and the dynamic mode decomposition (DMD) is discussed in Section 3.

In the present paper, we apply fully spatio-temporal Proper Orthogonal Decomposition to the turbulent channel flow. Wall turbulence is characterized by a variety of spatiotemporal scales interacting in a highly complex fashion [START_REF] Robinson | Coherent motions in the turbulent boundary layer[END_REF][START_REF] Dennis | Coherent structures in wall-bounded turbulence[END_REF]. 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 λ z+ ∼ 100 -150 in the spanwise direction and a length of 600 -1000 wall units [START_REF] Kim | The production of turbulence near a smooth wall in a turbulent boundary layer[END_REF], [START_REF] Stanislas | Vortical structures in the turbulent boundary layer: a possible route to a universal representation[END_REF], [START_REF] Jiménez | Near-wall turbulence[END_REF]), which are units based on the fluid viscosity and wallfriction velocity u τ , and will be denoted with a + throughout the paper. These streaks are associated with longitudinal vortical structures pushing low-speed fluid upwards and bringing high-speed fluid downwards, which results in a strong contribution to the Reynolds stress [START_REF] Kim | The production of turbulence near a smooth wall in a turbulent boundary layer[END_REF]. This contribution is highly intermittent in space and time, with 'bursts' of turbulence production. The time scale of the bursts is difficult to assess. Eulerian measures of the bursts yield time scales on the order of 300 wall units [START_REF] Blackwelder | Scaling of the bursting frequency in turbulent boundary layers[END_REF]. This value is based on the VITA (Variable Interval Time Average) criterion applied to the streamwise velocity and depends on the choice of a particular threshold, which makes it difficult to provide an absolute value for the bursting period. However the value appears to be independent of the Reynolds number. Minimal flow units, i.e. periodic domains of small extent [START_REF] Jiménez | Near-wall turbulence[END_REF]) were characterized by global characteristic time scales of 300-400 wall units, which appear to correspond to regeneration cycles of coherent structures [START_REF] Hamilton | Regeneration mechanisms of near-wall turbulence structures[END_REF], but relating this time scale with that of a full-scale turbulent flow is difficult.

Classic spatial POD has been previously used to investigate the boundary layer (e.g., see [START_REF] Aubry | The dynamics of coherent structures in the wall region of the wall boundary layer[END_REF], [START_REF] Moin | Characteristic-eddy decomposition of turbulence in a channel[END_REF], [START_REF] Podvin | A low-dimensional approach for the minimal flow unit[END_REF]). Spatial eigenfunctions are determined from second-order 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 spatio-temporal patterns.

In this paper we focus on a relatively moderate Reynolds number Re τ = 590, based on the fluid viscosity, channel half-height and friction velocity. At these moderate Reynolds numbers our focus will be on the wall layer as large scales such as those observed by [START_REF] Smits | High-reynolds number wall turbulence[END_REF] 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 spatio-temporal POD. Section 3 discusses the connections of our approach with snapshot POD and DMD. Section 4 presents the POD results while the contributions to the turbulent kinetic energy equation associated with each mode are examined in Section 5, followed by a conclusion in Section 6 which summarizes the key observations and results.

Description of the procedure

Full spatio-temporal POD

Owing to the homogeneity of the statistics in the spatial directions x and z and their stationarity, POD modes are Fourier modes in the horizontal directions as well as in time, [START_REF] Berkooz | The proper orthogonal decomposition in the analysis of turbulent flows[END_REF]. 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 pth component of the velocity field (u p (x, y, z, t)) in the physical space is then written as

u p (x, y, z, t) = l k f u p lkf (y)e 2πi(lx/Lx+kz/Lz+f t/T ) , (2.1)
where u p lkf (y) denotes the pth velocity component in the Fourier space corresponding to the streamwise wavenumber l, spanwise wavenumber k and frequency f . Proper Orthogonal Decomposition is then applied in the wall-normal direction for each triad (l, k, f ):

u p lkf (y) = n a n lkf φ n,p lkf (y), (2.2)
with, for any two 4-tuples (l, k, f, n) and (l, k, f, m), a n lkf a m * lkf = δ nm λ n lkf .

(2.3)

Here, • represents the usual ensemble average over the space of all possible flow realizations, δ mn is the Kronecker-delta function, n denotes a specific POD mode for each triad, and * denotes complex conjugation. Note that a n lkf is independent of p as the three velocity components are stacked into a single vector. The complex, stochastic coefficients a n lkf are uncorrelated and their variance is equal to λ n lkf . The eigenfunctions φ n,p lkf and eigenvalues λ n lkf are obtained by solving the following eigenproblem, where the autocorrelation is estimated from taking an ensemble average:

Y 0 u lkf (y)u * lkf (y ) φ n lkf (y )dy = λ n lkf φ n lkf (y).
(2.4)

Here, u lkf (y) represents the single velocity vector with the three components stacked in it. Upon discretization in y, φ n lkf (y) represents the single eigenvector containing the three components. By construction, the eigenvectors are orthonormal with respect to the Euclidean inner product, and λ n lkf can be interpreted as the energy content in each mode. Similarly, for each triad (l, k, f ), the modes (i.e. different n's) are sorted by energy. Since the left-hand-side of Equation (2.4) is a Hilbert-Schmidt integral operator, its eigenvalues are real and non-negative. The eigenvalues λ n lkf associated with the triplets (l, k, f ) can then be gathered and sorted globally according to their magnitude. We will call mode number the global index N (l, k, f, n) associated with the sorted modes over spatial wavenumbers (l, k), frequency (f ) and quantum number n, with lower mode number indicating a larger eigenvalue. In all that follows, we focus only on the most energetic modes.

Numerical implementation

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 T , each of which constitutes a sample. The samples are therefore not independent realizations, but we assume that the period T 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 T cannot be too large in order to allow for a reasonable number of samples N s to be constituted. We note that as in [START_REF] Moin | Characteristic-eddy decomposition of turbulence in a channel[END_REF] and [START_REF] Podvin | On the adequacy of the ten-dimensional model for the wall layer[END_REF], 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 n t velocity fields uniformly sampled at a rate δt. The value of δt is chosen such that the Nyquist frequency is well above the range of characteristic frequencies of the flow. It is however limited by computational tractability, as large values of n t will increase memory requirements for the samples. For most of the test cases, we have fixed the value of T to be 3000 wall units.

We use the finite-volume code SUNFLUIDH to simulate the incompressible flow in a channel. Further details on the configuration and numerical scheme can be found in [START_REF] Podvin | POD-based wall boundary conditions for the numerical simulation of turbulent channel flows[END_REF]. Periodic boundary conditions are enforced in the streamwise (x) and the spanwise direction (z), and the channel lengths along these directions are L x = 2πh and L z = πh, where h is the channel half-height. The total domain is discretized using a grid of size 256 × 256 × 256, with a uniform grid for the horizontal directions x and z and a hyperbolic tangent-distributed mesh for the wall-normal direction y. The velocity components in the streamwise, wall-normal and spanwise direction will be denoted respectively by u, v, w or equivalently by u 1 , u 2 , u 3 . The simulation is conducted at Re τ = 590, based on the friction velocity u τ and channel half-height h, which given block, N s denotes the number of samples, δt + is the time gap between successive snapshots, and T + is the total time period for a given block (expressed in viscous units).

Figure 1: Schematic of the methodology. Each grey slab on the left represents an instantaneous realization of the flow in the physical domain (the wall-normal direction corresponding to the thickness of the slab), while each slab on the right corresponds to the Fourier transform of a block of n t realizations in time as well as in the horizontal directions x and z. For each Fourier triad (l, k, f ) the autocorrelation tensor is computed by averaging over the N s blocks and an independent eigenproblem is solved.

corresponds to the simulation of [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Reτ = 590[END_REF]. The database consists of 16500 flow realizations, separated by δt + = 5.9. The time spanned by the database represents about 10 5 viscous time units.

In most of the paper, we focus on the domain y + < Y + = 80, but we also considered the full boundary layer where Y + = 590. 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 600-1000 wall units [START_REF] Jiménez | Near-wall turbulence[END_REF])), the streamwise extent of the domain was limited to L x = π/2h representing 900 wall units. The full spanwise extent of the domain (about 1800 wall units) was considered.

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 the configurations D1-D3. These datasets correspond to the wall layer, Y + = 80. The fourth one (D4) corresponds to the full boundary-layer height, Y + = 590.

Connections with alternative decomposition methods

We now briefly explore the connection of the spatio-temporal POD with the snapshot POD and the Dynamic Mode Decomposition, two established methods for representing and analyzing fluid flows. The numerical cost and requirements of the present method are also discussed.

Snapshot POD

Snapshot POD [START_REF] Sirovich | Turbulence and the dynamics of coherent structures[END_REF] relies on an ergodicity assumption to identify the ensemble average operator • with the time-average. In practice, one has to rely on a finite set of samples (snapshots) and the average operator reduces to an empirical algebraic average. Snapshot POD then considers the empirical spatial autocorrelation tensor K U (x, x ) = n -1 t nt j=1 u (x, t j ) u (x , t j ), with x the vector-valued continuous space variable. Autocorrelation tensor eigenfunctions φ j (x), associated with the largest eigenvalues λ j , are the dominant POD modes. The operator involved in the Fredholm equation of the snapshot POD method being Hilbert-Schmidt, its eigenfunctions are orthonormal with respect to the retained inner product.

The spatial autocorrelation tensor K U is estimated from a finite set of samples and then only provides an approximation of the true autocorrelation tensor

K U (x, x ) = lim |Ωt|→∞ |Ω t | -1 Ωt u (x, t) u (x , t) dt.
In general, the spatial POD modes are then only an approximation of the true autocorrelation eigenvectors. In the spatio-temporal POD approach, we take advantage of the closed-form solution of the Fredholm equation in homogeneous dimensions. No approximation is then introduced and the spatio-temporal modes are the exact eigenvectors of the autocorrelation tensor along the homogeneous dimensions.

From a dynamical perspective, the temporal evolution of snapshot POD mode j is given by the orthogonal projection of the flow field onto the corresponding eigenfunction: a j (t) = (φ j , u (t)), where u (t) ≡ u (x, t) and (, ) represents the (typically weighted Euclidean) inner product. The time evolution of mode j then inherits many properties of the flow field such as its wide frequency content. It is important to note that, since it results from a projection, there is no guarantee that the dynamics of a snapshot POD mode is smooth in time.

In contrast, the present spatio-temporal POD relies on a spectral decomposition in the homogeneous dimensions, including time. Each mode then follows a smooth (harmonic) dynamics. The original flow field being approximated by a finite linear combination of smooth modes, it remains smooth in time. More generally, spatio-temporal POD enjoys structure in both time and homogeneous dimensions and allows to interpret modes as coherent structures, localized in the frequency-spatial wavenumber space. In contrast, snapshot POD only enjoys structure in space and lacks structure in time, preventing identification of its modes with coherent structures which dynamics is essentially smooth.

DMD

The Dynamic Mode Decomposition (DMD) is another popular method for modal decomposition of fluid flows. Different definitions have been considered in the literature but we here focus on the original formulation discussed in [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. Considering the flow field at a collection of N x points in space as a vector-valued observable, and assuming it is a state vector for the underlying physical system, the DMD can be closely related to the Koopman theory, [START_REF] Mezić | Spectral properties of dynamical systems, model reduction and decompositions[END_REF]; [START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF]. In a nutshell, the DMD estimates the eigenvectors and associated eigenvalues of the linear operator mapping the discrete flow field at a given time to a subsequent time ∆t in the future. Specifically, n t time-ordered snapshots of the flow field sampled every ∆t in time are collected in a matrix U . Each snapshot is of size N x and constitutes a vector-valued observable u of the state vector of the underlying physical system. Letting U 1 be the first n t -1 columns of U and U 2 the last n t -1 columns, the DMD is concerned with the characterization of the linear operator A such that U 2 A U 1 . Formulating the approximation as an optimization problem, the eigendecomposition of the matrix A can be obtained: A Ψ = Ψ Λ, with Ψ ∈ C (nt-1)×(nt-1) the matrix of eigenvectors and Λ ∈ C (nt-1)×(nt-1) the diagonal matrix of eigenvalues λ j . The flow field u (x, t) is identified with u (t) and can then be approximated at a time t + l ∆t as

u (t + l ∆t) ≈ Ψ diag Ψ -1 u (t) λ l = Ξ λ l , ∀ l ∈ Z,
where

λ l = λ l 1 λ l 2 . . . λ l nt-1
T and the matrix Ξ = ξ 1 ξ 2 . . . ξ nt-1 contains the so-called DMD modes ξ j .

Assuming the collection of snapshots is linearly independent and subtracting the empirical algebraic time-average from the snapshot data U , the DMD is equivalent to the temporal Discrete Fourier Transform (DFT), [START_REF] Chen | Variants of dynamic mode decomposition: boundary condition, koopman, and fourier analyses[END_REF]. Specifically, denoting zero-mean quantities with a superscript , the vector-valued observable u (t + l ∆t) then obeys

u (t + l ∆t) = nt-1 j=1 exp 2πilj n t ξ j .
In this context, DMD decomposes the flow field in monochromatic spatial modes ξ j oscillating at a given frequency j/ (∆t n t ) in time. Both the spatio-temporal POD and DMD then derive modes oscillating at a given frequency. The sampling frequency being a multiple integer of these frequencies, both methods represent the data with modes at the same frequencies.

However, an essential difference between DMD and POD is that there is no notion of ensemble average in DMD (only one sample, corresponding to the set of snapshots, is considered), and DMD modes are extracted from mapping one snapshot in the time series to the next one. Standard DMD is applied to one flow realization, and thus a mode computed at given frequency would have statistical variability over a set of realizations. In contrast, in the present decomposition, ensemble average is a key feature of the procedure and requires several samples. POD modes are extracted from the autocorrelation tensor built from the different samples. They constitute an optimal basis that account for the variability of modes over an ensemble of realizations. We point out that in our numerical implementation, the different samples are obtained by breaking down a single series of snapshots into different non-overlapping blocks.

Requirements and cost of the decomposition

The spatio-temporal POD discussed in this paper relies on the same information as the snapshot POD or the DMD. Owing to the DFT in time and homogeneous spatial directions, it however involves distinctive differences such as the requirement of collecting time-resolved snapshots. This is in contrast with the other methods which only require time-sampling in pairs (DMD) or a set of snapshots sampled in time according to an approximately ergodic measure (snapshot POD).

Obtaining time-resolved data in the full spatial domain from a numerical simulation is not a problem, but requires integration over a sufficiently long period of time, especially if it needs to be broken down into different samples. In contrast, in an experiment, obtaining different samples over relatively long periods of time is relatively easy, but access to full spatial resolution may be more difficult, especially for time-resolved data, despite recent advances in PIV techniques (see for instance [START_REF] Westerwheel | Particle image velocimetry for complex and turbulent flows[END_REF]). However, it should be borne in mind that ensemble average and Fourier transform commute, so that full information in the spatial domain does not have to be acquired simultaneously, since the autocorrelation tensor for given spatial separations in homogeneous directions can be computed independently. If there is one homogeneous direction (say x), the autocorrelation tensor K U (x 1 , y 1 , x 2 , y 2 ) = f (x 2 -x 1 , y 1 , y 2 ) can be evaluated independently for each separation x 2 -x 1 from two simultaneous planes x = x 1 and x = x 2 : the information does not have to be acquired for all positions x at the same time, only pairs of simultaneous positions are necessary. For the same reason, in the case of two homogeneous directions (x and z), one can show that the full spatial autocorrelation tensor at all separations (x 1 -x 2 , z 1 -z 2 ) can be entirely recovered from two orthogonal planes x = cst and z = cst (as was done in the WALLTURB experiment, [START_REF] Stanislas | Progress in Wall Turbulence: Understanding and Modeling[END_REF].

We now briefly discuss the numerical cost of the solution method presented in this paper. We consider the common situation where the sampling is uniform in time and sufficient finely resolved with respect to the finest scales one is interested in. A Discrete Fourier Transform is then applied in homogeneous dimensions. We use a multidimensional Fast Fourier Transform (FFT) which essentially amounts to the composition of a sequence of one-dimensional FFTs along each homogeneous dimension. The numerical cost of this step then retains the O (N log (N )) scaling of FFT techniques, where N is here N l , N k or n t . Efficient algorithms also exist in case of non-uniform sampling, e.g., Ruiz-Antolín & Townsend (2018), and retain the O (N log (N )) scaling. The remaining step involves a POD in the non-homogeneous dimension for each atom of triad (l, k, f ). In the present case, it then reduces to a set of independent one-dimensional PODs.

Finally, we would like to stress that the eigenvalue problems involved in every method discussed here (spatio-temporal POD, spatial POD, DMD) can be advantageously solved by recent numerical techniques. In particular, sketching and randomized methods [START_REF] Halko | Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions[END_REF][START_REF] Kannan | Randomized algorithms in numerical linear algebra[END_REF] can very significantly alleviate the numerical cost and the memory requirement in cases only part of the spectrum is required, as is often the case in our applications.

Results from spatio-temporal POD

POD eigenvalues

Figure 2a 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 3n t × N l × N k × N y , where n t is the number of instantaneous fields contained in a sample, N l and N k are respectively the numbers of streamwise and spanwise wavenumbers, and N y is the number of grid points in the wall-normal direction. For D3 this corresponds to about 4 × 10 8 degrees of freedom. The largest eigenvalue N (l, k, f, n) = 1 where l = k = f = 0 and n = 1 corresponds to that of the mean mode (discussed in the next section). For each configuration, we have applied spatio-temporal POD to the data without removing the sample mean and compared the results to the case where the sample mean was removed. In both cases the eigenvalue distribution was very similar, and the deviations in eigenfunctions were negligibly small at least for the most energetic modes. In the rest of the paper we will focus on the n = 1 eigenfunctions. We note that the eigenfunctions corresponding to n = 1 capture about 80% of the turbulent kinetic energy, which is in agreement with [START_REF] Moin | Characteristic-eddy decomposition of turbulence in a channel[END_REF]'s results.

Figure 2a shows that the spectrum eigenvalues 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 eigenvalues in D1 and D2 are similar to that in D3 at low values of N , but have higher energy levels at large values of N , 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 λ(N ) ∼ N -1.27 . Unlike the other configurations, D4 corresponds to the full boundary layer 0 < y + < 590. 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 0 < y + < 80, which suggests self-similarity (see next section). Figure 2b shows the fraction of turbulent kinetic energy (TKE) captured by the first N modes for each N (as we are considering fluctuations, the contribution from the mean mode was set to zero). The fraction of TKE captured by the most energetic 100 mode numbers is relatively high, but increases only slightly for mode numbers higher than 1000, which highlights the complexity of the flow. The slow convergence in Fourier space also indicates that the structures are highly localized in space and time, owing to the fundamentally intermittent nature of turbulence [START_REF] Frisch | Turbulence[END_REF]. This is confirmed by Table 2, which shows that the first 200 eigenvalues capture about 21% of the TKE which is increased only to 32% when 5000 modes are included, for the D3 configuration. The trend remains the same for the D4 configuration, where 200 eigenvalues capture around 24% of TKE which increases only to 34% when 5000 eigenvalues are included.

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 ( l, k, f ), expressed in wall units, are related to Figure 3a shows how the sum of the eigenvalues over streamwise wavenumber and frequencies ( l f λ 1 lkf ) varies as a function of k+ . Results show a good agreement with the rescaled standard energy spectrum obtained by [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Reτ = 590[END_REF] at y + = 80. 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 n = 1 eigenvalues over all (l, k) modes (i.e. l k λ 1 lkf ) as a function of frequency f+ for the different configurations D1-D4. As observed previously, aliasing effects can be observed for the fields sampled in time at a lower rate (D1-D2), but the trends are very similar. A marked increase in the energy is observed in the frequency range of 0.003 -0.005, which corresponds to time scales of 200-300 viscous units. This value is in good agreement with the duration of the regeneration cycle identified by [START_REF] Hamilton | Regeneration mechanisms of near-wall turbulence structures[END_REF] and [START_REF] Jimenez | The minimal flow unit in near-wall turbulence[END_REF] in minimal domains, as well as the investigations of larger-scale domains reported in [START_REF] Jiménez | Near-wall turbulence[END_REF]. 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 f+ = 0.005, while it is located at f+ = 0.0044, in the region y + < 80.

As far as we know this is the first evidence of an objective time scale identified in the channel configuration. The time scale did not appear to be affected by spatial variability: it did not depend on the spanwise extent (different domain widths were considered) or the streamwise extent of the domain (two domains of different streamwise sizes were compared). We note that the dominant frequency of the nonlinear terms (e.g. the Reynolds stresses) should be on the order of 100-150 wall units. An interesting connection, which will need further exploration, can be made with wall reduction control schemes based on wall spanwise oscillations such as those of [START_REF] Quadrizio | Critical assessment of drag reduction through spanwise wall oscillations[END_REF] and [START_REF] Choi | Turbulent boundary layer control by means of spanwise wall oscillations[END_REF], where it was shown that the optimal oscillation period is about 100-150 wall units.

The predominance of the characteristic frequency fc+ at f+ = 0.0044 is confirmed by Table 3, which shows the top 30 eigenvalues (denoted by λ n lkf ) with the corresponding N l k f n λ n lkf 1 0 0 0 1 24.98 2 0 2 0 1 0.008133 3 0 1 0 1 0.004205 4 0 3 0 1 0.003754 5 0 4 0 1 0.001489 6 0 3 1 1 0.001277 7 0 2 1 1 0.001222 8 0 1 1 1 0.001128 9 0 4 1 1 0.000923 10 0 5 0 1 0.000655 11 0 3 13 1 0.000591 12 0 2 14 1 0.000552 13 0 5 0 1 0.000523 14 0 2 13 1 0.000505 15 0 3 12 1 0.000457 N l k f n λ n lkf 16 0 3 2 1 0.0004568 17 0 2 2 1 0.0004559 18 0 6 0 1 0.000451 19 0 4 12 1 0.000415 20 0 3 14 1 0.000412 21 0 4 2 1 0.000396 22 0 6 1 1 0.0003806 23 0 0 1 1 0.0003643 24 0 5 2 1 0.000349 25 0 4 11 1 0.0003475 26 0 2 12 1 0.0003445 27 0 4 13 1 0.000335 28 0 5 12 1 0.000314 29 0 2 15 1 0.000295 30 0 1 14 1 0.000289

Table 3: Top 30 most energetic eigenvalues λ n lkf along with their corresponding (l, k, f, n). Fourier mode indices (l, k, f ) in space and time and quantum mode n for the configuration D3. Although a large fraction (18 modes) of the first 30 modes are associated to f 2, which corresponds to long time scales that are outside the scope of our analysis, all the other modes are characterized by frequencies between 11 and 15, which correspond to the previously identified time scale of about 200-300 viscous units. The most energetic of these modes corresponds to f = 13, i.e. fc+ . We note that all the modes in the table are characterized by n = 1, which corresponds to the largest eigenvalue for a triad (l, k, f ), and l = 0, which corresponds to the streamwise-averaged flow.

To provide a comparison with standard analysis tools, Figures 4a andb show the temporal Fourier transform of the streamwise velocity component averaged along horizontal directions for two different heights: y + = 3 and 15. Although there seems to be some energy increase around the characteristic frequency fc+ , evidence of a local peak is not clear in these figures. This illustrates the usefulness of the new POD implementation.

Figures 5a and5b respectively show the distribution of the top 5000 eigenvalues along k+ and f+ 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 f+ at around 0.009 and 0.0135, which correspond to the harmonics of the frequency 0.0044.

Figure 6a shows that the energy spectrum integrated in frequency space f λ 1 lkf slowly decreases in an apparently self-similar manner with respect to both streamwise and spanwise wavenumbers. Figure 6b shows the corresponding variations in the frequency space of the energy spectrum integrated in the spanwise wavenumber space for different streamwise wavenumbers l. For the l = 0 mode, owing to spanwise reflection invariance, the plot is symmetric in the frequency space, and the characteristic frequency fc+ and its harmonics can be clearly identified. These peaks observed at l = 0 are still present in the l = 0 spectra, which confirms that the characteristic frequency f c+ is not an artifact of the streamwise average. However non-zero 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 [START_REF] Krogstad | Convection velocities in a turbulent boundary layer[END_REF]; [START_REF] Alamo | Estimation of turbulent convection velocities and corrections to taylor's approximation[END_REF]), but the question can be more easily addressed in the present framework where each POD eigenfunction is naturally associated with a phase velocity -f / l. The spectra k λ 1 lkf can be used to define a global convection velocity c at each streamwise wavenumber using

c( l) = - fmax l , (4.2)
where fmax = arg max

f k λ 1 lkf . (4.3)
Peaks in the spectra are located at f+ = (0.013, 0.027, 0.044) for respective streamwise wavenumbers of l+ = (0.0011, 0.0022, 0.0033). This corresponds to a global convection velocity of 12u τ , with a slight upward shift observed with increasing streamwise wavenumbers. This agrees well with the classical results [START_REF] Kreplin | Propagation of perturbations in the viscous sublayer and adjacent wall region[END_REF], [START_REF] Krogstad | Convection velocities in a turbulent boundary layer[END_REF], [START_REF] Wallace | Space-time correlations in turbulent flow: A review[END_REF], [START_REF] Alamo | Estimation of turbulent convection velocities and corrections to taylor's approximation[END_REF]).

Having access to the four-dimensional space provides a way to test the Taylor's frozen turbulence hypothesis [START_REF] Lumley | On the interpretation of temporal spectra in high intensity shear flows[END_REF], which states that the spatial spectrum in the streamwise direction E kx of u(x)u(x ) (where the overbar represents a spatial average) can be recovered from the temporal spectrum E f of u(t)u(t ) (the double overbar denotes here a temporal average) where f = kU and U is a suitable convection velocity. In the context of the present decomposition, this leads us to compare, for each spanwise wavenumber k z and quantum mode n, the temporal spectrum

Ẽ(f , k z , n) = l λ n lkf δk x ,
where f = f /T and k z = k/L z , with the streamwise spectrum for the streamwise wavenumber

k x E(k x , k z , n) = f λ n lkf δf ,
where k x = l/L x , using a suitable rescaling factor δk x /δf which can be obtained from

E = l k f λ 1 lkf δk x δk z δf = l k E(k x , k z , n)δk x δk z = f k Ẽ(f, k z , n)δf δk z .
For the first quantum mode n = 1, Figure 7 compares the spectra in the temporal direction E(k x , k z , 1) with the equivalent spectrum in the streamwise direction Ẽ(f, k z , 1). The frequency expressed in wall units f + is related to the the streamwise wavenumber k x+ such that k x+ U c+ = f + where we have taken U c+ = 12. For the range of corresponding frequencies, a good agreement can be observed between the spectra at each spanwise wavenumber k z .

POD eigenfunctions

Figure 8 shows that the reconstructed mean streamwise profile coincides with previous results [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Reτ = 590[END_REF]) and with the dominant mode (l = 0, k = 0, f = 0, n = 1) of the decomposition. Figure 9 shows the absolute value of each component of the most energetic mode φ p,n lkf , for ( l+ , k+ , f+ ) = (0, 0.0016, 0.0044), obtained for the datasets D1-D3. As expected, there is a good agreement between the different configurations, which Figure 8: Comparison of the mean streamwise profile reconstructed using the first eigenfunction for the mode pair (0, 0, 0) with that of [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Reτ = 590[END_REF], and with the mean profile obtained from the DNS.

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 10a,c and e show the shape of the most energetic eigenfunction components φ 1,p 0kf associated with the characteristic frequency fc+ ∼ 0.0044 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 wall-attached structures (Alamo et al. (2006), [START_REF] Podvin | On self-similarity in the inner wall layer of a turbulent channel flow[END_REF]). The values of the maxima are similar for both cross-stream components, also in agreement with previous observations [START_REF] Alamo | Self-similar vortex clusters in the turbulent logarithmic region[END_REF], [START_REF] Podvin | A few thoughts on proper orthogonal decomposition in turbulence[END_REF]), but the value tends to increase with the spanwise wavenumber for the streamwise and the wall-normal component, while it slightly decreases for the spanwise component. The monotonous evolution of the shape of the modes with the spanwise wavenumber suggests self-similarity, as was proposed in [START_REF] Podvin | A few thoughts on proper orthogonal decomposition in turbulence[END_REF].

Figures 10b,d and f compare the eigenfunctions obtained on the domain 0 < y + < 80 (for D3) and 0 < y + < 590 (for D4) for several spanwise wavenumbers. The eigenfunctions of D4 have been rescaled to have the same energy content as those of D3 on 0 < y + < 80. 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 wall-normal 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 [START_REF] Podvin | On self-similarity in the inner wall layer of a turbulent channel flow[END_REF]). Figure 11 shows the hierarchical organization of the eigenfunctions in various dimensions of the four-dimensional space. For each streamwise wavenumber and quantum mode, we define the intensity of the eigenfunctions as a function of height and absolute frequency (for the sake of clarity positive and negative frequencies are aggregated, although as seen in a previous section there is no symmetry):

| y + k = 1 k = 2 k = 3 k = 4 k = 5 k = 6 k = 7 k = 8 (
|φ 1,1 lkf | y + k = 1¸ ¿ k = 1¸ k = 3¸ ¿ k = 3¸ k = 5¸ ¿ k = 5¸ k = 7¸ ¿ k = 7¸ (
I ln (y, |f |) = k (λ n lkf |φ n lkf (y)| 2 + λ n lk-f |φ n lk-f (y)| 2 ).
It is represented in Figure 11 for the first quantum numbers n 4 and selected streamwise wanumbers l. Each quantum number n is characterized by n peaks at a given frequency. The wall-normal location of the highest peak increases nearly linearly with the quantum number n, and is about 10n in wall units. At a given wavenumber in the horizontal space, energy is therefore transferred from the more energetic to the less energetic modes towards the core region away from the wall. The range of frequencies associated with a high intensity increases only slightly with the quantum number n, but increases significantly with the streamwise wavenumber l with a shift towards higher frequencies, which is due to convection effects.

Figure 12 shows the flow generated by selected modes in figure 10 associated with fc+ for two different values of k+ . As these dominant modes have no streamwise variation (i.e. l = 0), the three-dimensional shapes of the modes are represented by a contour plot of the streamwise velocity component (u) and a velocity plot of the in-plane quantities (v, w). The modes are characterized by updrafts of low-speed fluid alternating with downdrafts of high-speed fluid, associated with vortical motions in the streamwise direction, in agreement with classic observations [START_REF] Corino | A visual investigation of the wall region in turbulent flow[END_REF]. The wallnormal extension of the structure decreases as the spanwise wavenumber increases, in agreement with Townsend's model [START_REF] Townsend | The Structure of turbulent shear flow[END_REF]) of wall-attached eddies (note the similarity of Figure 12 with Figure 9 in Jiménez ( 2013)). Figure 13 shows that the most energetic streamwise mode convected at a velocity of 12u τ also corresponds to high and low-speed streaks alternating in the streamwise as well as in the spanwise direction and associated with vortical motions in the cross-stream plane. The wall-normal position of the vortex centers depends on the vertical extension of the streaks. The contribution of the M most energetic modes to the fluctuating stresses can be evaluated using (4.4) where the subscript 'M ' is used to distinguish it as a reconstructed quantity using M most energetic modes. The bar indicates average over space and time and . denotes the usual ensemble operator. Note that the mean modes corresponding to (0, 0, 0, n) have been (1999). A non-negligible fraction of the turbulent intensity is recovered with the most energetic 200 modes: 60% for the streamwise velocity, 30% for the cross-stream 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 wall-normal component), which highlights the complexity of the flow.

u i u j M = N (l,k,f,n) M λ n lkf φ n, * i lkf φ m,j lkf ,

Energy budget of modes

It is of interest to determine the contributions made by the most energetic POD modes to the momentum equations. In pioneering studies, [START_REF] Hong | Application of large-eddy interaction model to channel flow[END_REF] and [START_REF] Gatski | Proper orthogonal decomposition based turbulence modeling[END_REF] examined the contributions of the dominant POD mode. Hong and Rubesin showed that the dominant fluctuating mode captured several essential characteristics of the fully turbulent flow. Gatski and Glauser considered the various contributions of the first POD eigenmode to the different terms of the momentum equation and showed that a turbulent model for the transport term could be built from one-mode estimates.

With the new variant of Proper Orthogonal Decomposition, projection of the Navier-Stokes 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 N (l, k, f, n).

For consistency with standard analysis, we decompose the pth component of the velocity field into its mean part U p which, as we have seen above, can be represented by the mode (0, 0, 0, 1) and a fluctuating part u p .

We rewrite the equations as

E p = - ∂u p ∂t -u j ∂U p ∂x j -U j ∂u p ∂x j -u j ∂u p ∂x j - ∂p ∂x p + ν ∂ 2 u p ∂x j ∂x j = 0, (5.1)
where U 1 = U (mean streamwise velocity) and U j = 0 if j = 2, 3 and ∂U p ∂xj = 0 only if p = 1 and j = 2.

To obtain a budget for each mode (l, k, f, n) we project Equation (5.1) onto the corresponding mode, i.e. we take the Fourier transform in the (l, k, f ) space, take the inner product in the wall-normal direction, and apply an ensemble average. This gives where * denotes complex conjugation.

• The first term -∂u p ∂t yields a contribution -I n lkf = -if λ n lkf . 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:

P n lkf = - Y 0 u 2 dU dy a * n lkf φ * n,1 lkf dy = -λ n lkf Y 0 dU dy φ n,2 lkf φ * n,1 lkf dy.
(5.

3)

The third term represents the convection effect of the mean field:

C n lkf = - Y 0 U ∂u p ∂x 1 a * n lkf φ * n,p lkf dy = - Y 0 m U ∂φ m,p lkf ∂x 1 φ * n,p lkf a m lkf a * n lkf dy = -i lλ n lkf Y 0 U φ n,p lkf φ * n,p lkf dy.
(5.4)

• The last term on the left-hand-side of the equation corresponds to the viscous diffusion term D n lkf defined as:

D n lkf = λ n lkf Y 0 ∂ 2 φ n,p lkf ∂x j ∂x j φ * n,p lkf dy = λ n lkf (-l2 -k2 ) Y 0 φ n,p lkf φ * n,p lkf dy + Y 0 d 2 φ n,p lkf
dy 2 φ * n,p lkf dy .

(5.5) Integrating the last term by parts, one has

D n lkf = (-l2 -k2 ) - Y 0 dφ n,p lkf dy dφ * n,p lkf dy dy + dφ n,i lkf (Y ) dy φ * n,i lkf (Y ).
(5.6)

Except for the last term corresponding to a boundary effect, all contributions to D n lkf are real and negative, which corresponds to an energy loss for the mode, as expected.

Note that I n lkf , P n lkf , C n lkf and D n lkf only require information about the mode (l, k, f, n) (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 a n lkf . • As in the classic derivation [START_REF] Aubry | The dynamics of coherent structures in the wall region of the wall boundary layer[END_REF]) , 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

P r n lkf = p lkf (Y )a * n lkf φ * n,2 lkf (Y )
, where p lkf (Y ) is the Fourier transform (along x, z and t) of pressure at height Y , 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 [START_REF] Aubry | The dynamics of coherent structures in the wall region of the wall boundary layer[END_REF].

• Finally, the quadratic interactions can be evaluated as

Q n lkf = - l k f m p a m l k f a p l-l k-k f -f a * n lkf Y 0 ∂φ m,i l k f ∂x j φ p,j l-l k-k f -f φ * n,i lkf dy.
(5.7) Q n lkf requires information about the triple correlations of the coefficients a n lkf , and involves first derivatives of the eigenfunctions. This term characterizes the energy transfer from the different modes to the mode (l, k, f, n). It can be seen as a non-isotropic extension of the energy transfer function defined in isotropic turbulence [START_REF] Zhou | Interacting scales and energy transfer in isotropic turbulence[END_REF], and involves triads of modes in the (l, k, f ) space. It also corresponds to the forcing term of the resolvent analysis [START_REF] Mckeon | The engine behind (wall) turbulence: perspectives on scale interactions[END_REF]. If all other terms are known, the quadratic interaction terms Q n lkf can be determined from the budget equation:

I n lkf = P n lkf + C n lkf + D n lkf + P r n lkf + Q n lkf .
(5.8)

The real part of the equation represents a balance between the production P n lkf and dissipation term D n lkf , which depend only on the characteristics of the mode (l, k, f, n), and the interaction terms due to pressure P r n lkf and convection by velocity fluctuations Q n lkf , 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 f of the mode (I n lkf ) with the different physical mechanisms. Figure 15 represents the contributions of the different terms to the equations for the largest 200 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 (l = 0) and a quantum number n = 1. All terms are evaluated directly, except the quadratic term Q 1 0kf which is evaluated using Equation (5.8). 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 higher-order contributions were significant, which is typical of the closure problem. We note that the x-axis (N) represents modes ordered by energy, so that continuity in wavenumber or frequency space is not enforced, which may explain the "noisy" appearance of the plots.

For the sake of clarity, the left plots show mode number N from 1 -30 and the right plots show N from 30 -200. The top row (Figures 15a andb) shows that the real part of the production term P n lkf is essentially balanced by the sum of the real part of the dissipation and the quadratic terms. For the less energetic modes N > 30 (Figure 15b) 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 15 shows that the pressure interaction term is very small compared to the other terms, in agreement with [START_REF] Aubry | The dynamics of coherent structures in the wall region of the wall boundary layer[END_REF]'s derivation where it is modeled as a stochastic term of small amplitude. If we neglect the influence of the pressure term we have

-Re[Q 1 0kf ] ∼ P 1 0kf + D 1 0kf .
(5.9)

Figures 15c andd show that for the most energetic modes, as l = 0, the frequency f of the mode is directly related to the imaginary part (denoted by Im) of the nonlinear contributions to the mode (l, k, f, n):

Im[Q 1 0kf ] ∼ i f λ n 0kf .
(5.10)

This corresponds to the following closure approximation: (5.11) 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 Q n lkf can be determined from Equation (5.8), the Q n lkf φ n,p lkf (y)e 2πi(lx/Lx+kz/Lz+f t/T ) (5.12) Substitution of the velocity decomposition into the left-hand side of the equation requires performing a cumbersome convolution on all coefficients a n lkf , but the right-hand-side provides a straightforward expression of the quadratic terms in the stochastic POD basis of coefficients a n lkf , where coordinates are solely determined from second-order statistics. At each scale level (l, k, f, n), the effect of the quadratic terms is to stretch and to rotate the corresponding velocity mode by a factor Q n lkf . Direct modeling of the distribution of a n lkf could therefore lead to new formulations of turbulence models in the decomposition framework.

-Q 1 0kf ∼ l k f m p a m l k f a p l-l k-k f -f a * n 0kf ∂φ m,i l k f ∂x j φ p,j l-l k-k f f φ * 1,i 0kf dy ∼ P 1 0kf + D 1 0kf -i f λ 1 0kf .

Conclusion

Spatio-temporal 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 straightforward 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, which could have significant implications for control. Convection velocities could be directly defined from the POD spectrum. A global convection velocity on the order of 12u τ was identified in the wall layer, in good agreement with previous approaches. Examination of the spectrum provided an assessment and validation of the Taylor's frozen turbulence hypothesis.

About 30% of the turbulent kinetic energy was captured by the 200 most energetic modes. The most energetic modes were found to have a self-similar shape that appeared largely independent from the wall-normal extent of the decomposition domain, which shows the coherence of the motions over the height of the boundary layer. The modes appeared to be hierarchically organized, with a number of peaks and a maximal peak location at l = 0 directly proportional to the quantum number n, indicating that energy cascading towards the higher-order modes is directed away from the wall into the core region.

Finally, substitution of the decomposition into the Navier-Stokes 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 spatio-temporal Proper Orthogonal Decomposition.

Figure 2

 2 Figure 2: a) Log-log plot of the top 5000 eigenvalues for D1-D4 configurations. An asymptotic fit for the decay rate of the spectra for D3 is also shown. b) Percentage of turbulent kinetic energy (TKE) captured by the N most energetic modes for D3.

Figure 3

 3 Figure 3: a) Integrated energy spectrum ( l f λ 1 lkf ) as a function of k+ ; the data is compared with Moser et al. (1999)'s energy spectrum at y + = 80 (with appropriate rescaling). b) Integrated energy spectrum ( l k λ 1 lkf ) as a function of f+ for 0 < y + < 80 (D1-D4).

Figure 4 :

 4 Figure 4: Horizontally-averaged (x and z) temporal Fourier transform of the streamwise velocity component u as a function of f+ in a log-log scale at a) y + = 3 and b) y + = 15. A vertical line is shown at f+ = 0.005 to indicate the peak in both figures.

Figure 5 :Figure 6

 56 Figure 5: Eigenvalue distribution of the top 5000 eigenvalues for D3 configuration as a function of a) f+ and b) k+ .

Figure 7

 7 Figure 7: left: Eigenvalue spectrum l λ 1 lkf δk x+ as a function of f + ; right: Eigenvalue spectrum f λ 1 lkf δk x+ as a function of k x+ U c+ , where U c+ = U = 12. The colorbar is based on the decadic log scale.

Figure 9 :

 9 Figure 9: Comparison of the eigenfunctions φ 1 lkf corresponding to l+ = 0, k+ = 0.0016, f+ = 0.0044 for different configurations.

Figure 10 :

 10 Figure 10: Dominant eigenfunctions φ 1 0kf corresponding to fc+ and for increasing values of k+ . Top row: Streamwise component; Middle row: Wall-normal component; Bottom Row: Spanwise component; Left column: Comparison for 0 < y + < 80 (D3) and different spanwise wavenumbers (k = 1 to 8 which corresponds to k+ between 0.0011 to 0.0086)the legend for a), c) and e) is given in a); Right column: Comparison between the domains 0 < y + < 80 (D3) and 0 < y + < 590 (D4) for selected wavenumbers k = 1, 3, 5, 7.

Figure 11 :

 11 Figure 11: Structure intensity I ln (y, |f |) = k λ n lkf |φ n lkf (y)| 2 + λ n lk-f |φ n lk-f (y)| 2 as a function of the vertical direction y and of the frequency |f | for the first quantum numbers n 4 and different streamwise wavenumbers l. From left to right: increasing quantum numbers n = 1, 2, 3, 4. From top to bottom: increasing streamwise wavenumbers l = 0, 2, 4, 6.

Figure 12 :

 12 Figure 12: Contour plot of streamwise velocity component reconstructed from eigenfunction in (z, y) plane along with a vector plot for (w, v) velocity components for the non-trivial dominant eigenmode corresponding to a) l = 0, k = 3, f = 13, n = 1 and b) l = 0, k = 16, f = 13, n = 1.

Figure 13 :

 13 Figure 13: Reconstruction of the most energetic mode corresponding to the peak in the frequency spectrum in figure 6b for l = 1, k = 16, f = -40, n = 1. Isosurfaces of streamwise component in (x, y, z) are shown along with two slice views at x/h = 0.15 and x/h = 0.35 showing a contour plot of the streamwise component and a vector plot of the spanwise and wall-normal components. Units are arbitrary.

Figure 14 :

 14 Figure 14: Reconstruction of the turbulent intensities and Reynolds stress using 200 and 3000 modes. Comparison with Moser et al. (1999)'s and our DNS data is also shown. The legend is indicated in figure 14 a).

E

  p a * n lkf φ * n,p lkf (y)e -2πi(lx/Lx+kz/Lz+f t/T ) dydtdzdx, (5.2)

Figure 15 :

 15 Figure 15: Energy budget for the first 200 most energetic spatially fluctuating modes.

Table 1 :

 1 Description of different test cases; n t denotes the number of snapshots in a

	Name Lx/h Lz/h Y /h nt Ns δt+ T+
	D1	π/2	π	0.14 100 60 29.5 2950
	D2	π/2	π	0.14 100 120 29.5 2950
	D3	π/2	π	0.14 500 132 5.9 2950
	D4	π/2	π	1 500 132 5.9 2950

Acknowledgements

We are grateful to W.K. George for motivating discussions. We thank the anonymous referees for pointing out valuable references and making helpful suggestions. This work was supported by the Center of Data Science from the Paris-Saclay University. Computations were carried out at IDRIS-GENCI (project 02262).