domain: spi phys

text:
Sean Symon
email: ssymon@unimelb.edu.au

Denis Sipp

Beverley J Mckeon

A tale of two airfoils: resolvent-based modelling of an oscillator vs. an amplifier from an experimental mean

Keywords: 1, Introduction Recent studies, e, g

The flows around a NACA 0018 airfoil at a chord-based Reynolds number of Re = 10250 and angles of attack of α = 0 • and α = 10 • are modelled using resolvent analysis and limited experimental measurements obtained from particle image velocimetry. The experimental mean velocity fields are data-assimilated so that they are solutions of the incompressible Reynolds-averaged Navier-Stokes equations forced by Reynolds stress terms which are derived from experimental data. Resolvent analysis of the dataassimilated mean velocity fields reveals low-rank behaviour only in the vicinity of the shedding frequency for α = 0 • and none of its harmonics. The resolvent operator for the α = 10 • case, on the other hand, identifies two linear mechanisms whose frequencies are a close match with those identified by spectral proper orthogonal decomposition. It is also shown that the second linear mechanism, corresponding to the Kelvin-Helmholtz instability in the shear layer, cannot be identified just by considering the time-averaged experimental measurements as an input for resolvent analysis due to missing data near the leading edge. For both cases, resolvent modes resemble those from spectral proper orthogonal decomposition when the resolvent operator is low-rank. The α = 0 • case is classified as an oscillator and its harmonics, where the resolvent operator is not lowrank, are modelled using parasitic modes as opposed to classical resolvent modes which are the most amplified. The α = 10 • case behaves more like an amplifier and its nonlinear forcing is far less structured. The two cases suggest that resolvent-based modelling can be achieved for more complex flows with limited experimental measurements.

 (2018)have been able to estimate the spectra at various points of the flow using only the aforementioned measurements. The success of the method relies on the low-rank nature of the (linear) resolvent operator, which one obtains after linearising the Navier-Stokes equations around the (turbulent) mean flow, and treating the nonlinear terms as a source of intrinsic forcing [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF]. The exact form of the forcing can often be neglected due to a dominant amplification mechanism at particular spatial wavenumbers and frequencies. Linear analysis alone does not reveal the amplitude and phase of the most amplified disturbances yet they can be calibrated using limited unsteady data (Gómez et al. 2016a;[START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF].

One of the major advantages of this approach is that time-resolved data at all locations in the flow field are not required. From the governing equations, it is possible to obtain a prediction of the temporal Fourier mode in the actual fluctuation field. The goal of Dynamic Mode Decomposition (DMD) [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF][START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF]) is also to compute modes which oscillate at a distinct frequency, suggesting that resolvent analysis is attempting to approximate DMD modes from the equations instead of time-resolved data. A striking similarity between the most amplified resolvent mode and the DMD mode was noted by [START_REF] Gómez | On the origin of frequency sparsity in direct numerical simulations of turbulent pipe flow[END_REF] in turbulent pipe flow and the relationship between resolvent analysis and DMD was formalised in [START_REF] Sharma | Correspondence between Koopman mode decomposition, resolvent mode decomposition, and invariant solutions of the Navier-Stokes equations[END_REF].

Connections between resolvent analysis and other data-driven approaches have been discussed by [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF]. In particular, they examined modes extracted by Spectral Proper Orthogonal Decomposition (SPOD), introduced by [START_REF] Lumley | Stochastic tools in turbulence[END_REF], which produces a basis for the fluctuation field at each temporal frequency and optimally represents the space-time flow statistics. [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF] determined that SPOD modes and resolvent modes are equivalent if the nonlinear forcing is white in space and time. In practice, this might not be the case since the resolvent framework is a closed loop [START_REF] Mckeon | Experimental manipulation of wall turbulence: A systems approach[END_REF][START_REF] Mckeon | The engine behind (wall) turbulence: perspectives on scale interactions[END_REF] where the nonlinear interactions between the velocity fluctuations constitute the forcing to the resolvent operator. Despite this fact, resolvent analysis is quite often able to predict the structure of energetic velocity fluctuations. Moreover, if the singular values of the resolvent operator are well separated, it often suggests a similar separation in the SPOD eigenspectra as illustrated by [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF] for a turbulent jet.

Resolvent analysis may be used for predicting structures even if the singular values are not well separated. Recent work by [START_REF] Rosenberg | The role of parasitic modes in nonlinear closure via the resolvent feedback loop[END_REF] demonstrated that considering a small subset of triadic interactions of 'correct' resolvent modes (i.e. modes computed at frequencies where the singular values are well separated) can approximate the nonlinear forcing, a process reminiscent of weakly nonlinear analysis [START_REF] Sipp | Global stability of base and mean flows: a general approach and its applications to cylinder and open cavity flows[END_REF] without being limited to the vicinity of a critical Reynolds number. When these interactions are fed back through the resolvent, they produce the correct structure. These modes were called parasitic in the sense that they feed off nonlinear interactions of other modes rather than emerge naturally from a linear amplification mechanism embedded in the operator.

The majority of the aforementioned studies have utilised a mean obtained from simulations although a growing number have considered experimental means [START_REF] Sasaki | Real-time modelling of wavepackets in turbulent jets[END_REF][START_REF] Beneddine | Unsteady flow dynamics reconstruction from mean flow and point sensors: an experimental study[END_REF]Symon & McKeon 2018;[START_REF] He | Data assimilation and resolvent analysis of turbulent flow behind a wall-proximity rib[END_REF]. The objective here is to obtain a reduced-order representation using a partially known mean obtained from experiments. We use data-assimilation to recover the mean velocity field on a larger domain and enforce no-slip boundary conditions on the airfoil surface. Data-assimilation is a process whereby experimental measurements are merged with computational fluid dynamics (CFD) (see [START_REF] Hayase | Numerical simulation of real-world flows[END_REF] and a review of recent work can be found in da [START_REF] Da Silva | An EnKF-based flow estimator for aerodynamic flows[END_REF]. The particular framework used here was developed by [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean flow reconstruction[END_REF] and later applied to flow around an idealised airfoil by [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF].

In this article, we consider the flow around a NACA 0018 airfoil at a chord-based Reynolds number of Re = 10250 at two angles of attack: α = 0 • (A0 case) and α = 10 • (A10 case). The flows are obtained experimentally from particle image velocimetry (PIV). As will be discussed later in the article, the A0 case is an oscillator-type flow while the A10 case shares characteristics of both an amplifier and oscillator. An oscillator has intrinsic dynamics which are insensitive to background noise while an amplifier filters and amplifies upstream noise in the downstream direction [START_REF] Huerre | Hydrodynamics and nonlinear instabilities[END_REF][START_REF] Sipp | Dynamics and control of instabilities in open-flows: a linearized approach[END_REF]. The resolvent norm, or largest singular value of the resolvent operator, tends to have sharp peaks isolated at resonant frequencies in the cases of oscillators, e.g. Symon et al. (2018). Amplifier flows, on the other hand, may have separation between the first and second singular values over a broad ranges of frequencies, e.g. [START_REF] Dergham | Stochastic dynamics and model reduction of amplifier flows: the backward facing step flow[END_REF]; [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]. Of particular interest to us is the nature of the nonlinear forcing for each type of flow. If the nonlinear forcing is structured, then it can be deduced from a limited number of nonlinear interactions between coherent structures, e.g. [START_REF] Rosenberg | The role of parasitic modes in nonlinear closure via the resolvent feedback loop[END_REF]. In the case of a turbulent jet [START_REF] Towne | Stochastic and nonlinear forcing of wavepackets in a Mach 0[END_REF], which is more representative of an amplifier, the nonlinear forcing field resembled incoherent turbulent fluctuations.

The principal objectives of this article can be summarised as follows:

(1) does dataassimilation of the mean velocity field make a measurable difference on predicting the velocity fluctuations? (2) what does SPOD reveal when applied to the nonlinear forcing?

(3) can additional structure be extracted from resolvent analysis if the singular values are not separated but the nonlinear forcing is approximated? ( 4) is there a fundamental difference between oscillator and amplifier type flows? (5) how far can we go in modelling with just a few mean flow measurements and some unsteady point measurements? These questions are addressed in the rest of the paper which is organised as follows. In §2, we describe the PIV experiments and in §3 we discuss the mathematical background of the analysis tools used to construct the resolvent-based model. In §4, the mean velocity fields are data-assimilated onto a larger, more resolved mesh. In §5, we apply SPOD to the velocity fluctuations and nonlinear forcing of both airfoils to highlight the low-rank behaviour in each flow. In §6, the data-assimilated mean velocity fields are used as an input to resolvent analysis and the results are compared to those of SPOD. We also try to establish new links between resolvent analysis and SPOD which depend on whether the flow is an oscillator or an amplifier. In §7, we use resolvent modes and, where appropriate, parasitic modes to construct a model of the coherent fluctuations in the two airfoil flows before concluding in §8.

Experimental setup

The flows around a NACA 0018 airfoil at a Reynolds number of Re = 10250 and two angles of attack, α = 0 • (A0) and α = 10 • (A10), are measured using PIV. The experiments are performed in a free-surface water facility with a test section measuring 1.6 m in length, 0.46 m in width and 0.5 m in height. The NACA 0018 airfoil is mounted vertically in the tunnel, so that its span is parallel to the test section height, and has a chord length of 10 cm. Its spanwise length is 48 cm, resulting in an aspect ratio of AR = 4.8. The laser sheet is provided by a YLF dual cavity solid-state laser and is centred at a height of 220 mm. [START_REF] Dunne | Dynamic stall on vertical axis wind turbine blades[END_REF] measured the out-of-plane velocity component to be less than 4% of the free-stream velocity, allowing the mean flow at this height to be considered two-dimensional.

The PIV setup consists of two Phantom Miro 320 cameras with 50 mm focal length Nikon lenses and 1:1.8 aperture. They have an overlap of 18% in the streamwise direction and sample the flow at a frequency of 125 Hz, corresponding to a time resolution 0.0082c/U . The camera resolution is 1920 × 1200 pixels and they are calibrated at 8.2 px/mm, resulting in a spatial resolution of 50 vectors per chord length or 0.02c. The seeding particles are hollow glass spheres (reference 110P8 with an average diameter of 11.7 µm and a specific gravity of 1.1) and the seeding density is approximately 0.1 particles per square pixel. To reduce the impact of surface reflections, the image intensity is calibrated using white-image subtraction and background-image subtraction, which are taken before each run and averaged over 100 snapshots. LaVision's software package DaVis computes the velocity vectors where a standard cross-correlation technique via Fast Fourier Transformation is applied to each sequential image. The window-size is reduced from 32 × 32 px 2 to 16 × 16 px 2 over three passes. Once the data are post-processed, missing vectors are interpolated using an average of all non-zero neighbourhood vectors and a median filter [START_REF] Westerweel | Universal outlier detection for PIV data[END_REF] is used for outlier detection. A sample of the post-processed data for both cases is displayed in figure 1. The mean velocity field, obtained after time-averaging two runs of 3,499 snapshots, is subtracted so that the coherent structures can be more readily identified. As seen in figure 1, velocity vectors below the airfoil cannot be obtained since the laser sheet is obstructed in this region. For the A0 case, the mean velocity field is symmetric so the mean flow can be reflected over the centreline (the transverse velocity changes sign) to obtain the mean below the centreline and in the shadow region. This is not possible for the A10 case, so the mean velocity fields for α = ±10 • are obtained separately and stitched together to obtain data around the entire airfoil.

Analysis tools

Data-assimilation of the mean velocity fields

The experimental means are data-assimilated using the framework devised by [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean flow reconstruction[END_REF] so that they can later be used as an input to resolvent analysis. Only the equations are reviewed here as its application to experimental data is discussed by [START_REF] Symon | Data assimilation of mean velocity from 2D PIV measurements of flow over an idealized airfoil[END_REF]. The flow is governed by the incompressible Navier-Stokes equations, which are non-dimensionalised by the free-stream velocity U ∞ and the airfoil chord c:

∂ t u + u • ∇u = -∇p + Re -1 ∇ 2 u, (3.1a) ∇ • u = 0, (3.1b)
where u and p are the velocity and pressure, respectively. We use a Cartesian coordinate system where u, v, and w are the streamwise velocity, transverse velocity, and spanwise velocity in the x, y, and z directions. Substituting u = u + u , where an overbar denotes a mean and a prime denotes a fluctuation, into (3.1) and averaging in time yields the incompressible Reynolds-averaged Navier-Stokes (RANS) equations:

u • ∇u + ∇p -Re -1 ∇ 2 u = f = -∇ • R, (3.2a) ∇ • u = 0. (3.2b)
The right-hand side of (3.2a), or f , represents the divergence of the Reynolds stress tensor R.

The data-assimilation algorithm treats f as an unknown and leverages partial knowledge of u to reconstruct not only f but also u on a larger domain which obeys prescribed boundary conditions. A constrained optimisation problem is formulated as follows

L(u † , p † , u, p, f ) = E(u) -u † , u • ∇u + ∇p -Re -1 ∇ 2 u -f -p † , ∇ • u , (3.3)
to minimise the difference between a measurement of u, which satisfies (3.2), and the experimental mean velocity field u exp measurement. In (3.3), (•) † signifies an adjoint quantity and a • b = Ω a • bdΩ. The scalar E quantifies the discrepancy and is given by

E(u) = 1 2 m -M(u) 2 M . (3.4)
The operator M projects the numerical data to the subspace spanned by the known PIV measurements m.

The variations of L with respect to each dependent variable are set to zero to obtain the equations which are solved iteratively until E is minimised. These include

-u • ∇u † + u † • ∇u T -∇p † -Re -1 ∇ 2 u † = δE/δu, (3.5a)
∇ • u † = 0, (3.5b) ∇ f E = u † . (3.5c)
and (3.2), which one obtains from variations with respect to the adjoint variables. Equations (3.5a) and (3.5b) are the adjoint incompressible RANS equations and (3.5c) is the steepest descent direction. It should be noted that only the solenoidal component of f is recovered by the algorithm (see [START_REF] Symon | Mean and unsteady flow reconstruction using data-assimilation and resolvent analysis[END_REF]) so the spanwise component of ∇ × f can be compared between the experiment and the assimilation. The direct and adjoint RANS equations are solved in FreeFem++ [START_REF] Hecht | New development in FreeFem++[END_REF]) on the computational domain Ω ∈ -5

x 15 ∪ -5 y 5 with the leading edge of the airfoil being centred at the origin. The equations are spatially discretised using quadratic basis functions for the velocity and linear basis functions for the pressure resulting in approximately 360,000 and 1,000,000 degrees of freedom for the A0 and A10 cases, respectively.

Spectral proper orthogonal decomposition

The velocity fluctuations obtained from the processed PIV datasets are analysed using two modal decomposition techniques. The first, introduced by [START_REF] Lumley | The structure of inhomogeneous turbulent flows[END_REF][START_REF] Lumley | Stochastic tools in turbulence[END_REF], is spectral proper orthogonal decomposition (SPOD). From the raw, experimental data, SPOD calculates energy-ranked modes which are orthogonal to one another and oscillate at a single frequency. The SPOD modes are computed using the procedure described in [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF] and the code developed by [START_REF] Schmidt | An efficient streaming algorithm for spectral proper orthogonal decomposition[END_REF]. A brief summary of the steps is presented below.

The PIV data are rearranged into the matrix

Q = q 1 q 2 • • • q n ∈ R m×n , (3.6)
where m represents the number of states, i.e. the streamwise and transverse velocity of all PIV vectors, and n is the total number of snapshots. Thus, the overall size of Q is 49,104 total states by 6,800 snapshots. Following Welch's method [START_REF] Welch | The use of fast Fourier transform for the estimation of power spectra: a method based on time averaging over short, modified periodograms[END_REF], Q is split into p smaller, overlapping segments, or blocks, containing n b snapshots where n b < n

Q (p) = q (p) 1 q (p) 2 • • • q (p) n b ∈ R m×n b . (3.7)
Each row of Q (p) is then Fourier-transformed in time using the discrete Fourier transform (DFT), yielding an ensemble of Fourier realisations

Q(p) = q(p) ω1 q(p) ω2 • • • q(p) ωn b ∈ C m×n b . (3.8)
The columns of Q(p) represent a Fourier mode instead of a snapshot in time as they did in Q (p) . The Fourier modes for a specific frequency ω can be arranged into a new data matrix Q(ω) = q(1)

ω q(2) ω • • • q(b) ω ∈ C m×b , (3.9)
where b is the number of blocks.

The cross-spectral density matrix for a specific temporal frequency Ŝ(ω) is

Ŝ(ω) = Q(ω) Q * (ω) ∈ R m×m . (3.10)
The SPOD eigenvectors V (ω) and eigenvalues Λ(ω) can be obtained via an eigenvalue decomposition of the cross-spectral density matrix

Ŝ(ω) V (ω) = V (ω)Λ(ω). (3.11)
Since ω is specified, the modes identified by SPOD are complex but do not come in pairs as they would using proper orthogonal decomposition (POD) or DMD, for example. The second mode in the pair can be calculated by performing SPOD for -ω.

In practice, the number of blocks is much smaller than the number of states in the problem, i.e. b m, so it is computationally less demanding to consider the following eigenvalue problem instead

Q * (ω) Q(ω) Θ(ω) = Θ(ω)Λ(ω), (3.12) as Q * (ω) Q(ω) ∈ R b×b ,
and compute the eigenvector of the full system using

V (ω) = Q(ω) Θ(ω)Λ(ω) -1/2 . (3.13)
It should be noted that the choice of ω is limited by n and the sampling frequency which, for both experiments, is 125Hz. As will be seen in the next section, this is a major disadvantage of SPOD, which requires a lot of data to obtain converged modes and high frequency resolution, in comparison to resolvent analysis. The PIV snapshots are divided into blocks containing n = 1700 snapshots with 50% overlap, resulting in b = 7 realisations of the flow. The resulting frequency resolution is ∆ω = 0.46. SPOD is performed using a rectangular windowing function (nearly identical results are obtained using a Hamming window). The choice of n and b is a compromise between suitable resolution of the lower frequencies, where the energy is most concentrated, and obtaining a sufficient number of realisations to identify low-rank behaviour. We find that the trends are unaffected as long as we consider at least 3 realisations of the flow.

Resolvent analysis

Once the mean flow has been data-assimilated, it can be used as an input to resolvent analysis, which has been shown to provide an efficient basis for velocity fluctuations by [START_REF] Mckeon | A critical-layer framework for turbulent pipe flow[END_REF] and many others. Subtracting (3.2) from (3.1) results in the equations governing the fluctuations

∂ t u + u • ∇u + u • ∇u + ∇p -Re -1 ∇ 2 u = -u • ∇u + u • ∇u = f , (3.14a) ∇ • u = 0. (3.14b)
The left-hand side of (3.14a) contains the linear terms and the right-hand side is lumped into a single forcing term f to account for the nonlinear terms. Substituting f = f e iωt , where ω is a real-valued frequency, implies u = ûe iωt and leads to the following inputoutput form of the linearised Navier-Stokes equations

û = C T (iωB -L) -1 C f = H(ω) f . (3.15)
The matrix H(ω) is the resolvent operator and the linear operators C, B, and L are defined as

C = 1 0 , B = 1 0 0 0 , L = -u • ∇() -() • ∇u + Re -1 ∇ 2 () -∇() ∇ • () 0 .
(3.16) L is the linear Navier-Stokes (LNS) operator, B restricts the nonlinear forcing to the velocity subspace, and C can be adjusted to restrict the spatial domain of the input and output vectors.

The singular value decomposition (SVD) of the resolvent operator yields basis functions for the nonlinear term and velocity fluctuations at each temporal frequency ω

H(ω) = ∞ j=1 Ψ (ω)Σ(ω)
Φ * (ω).

(3.17)

The matrices Φ(ω) and Ψ (ω) contain the optimal forcing and response modes, respectively, which are ranked by their kinetic energy gains contained in the diagonal matrix Σ(ω). They can also be interpreted as bases for the nonlinear forcing and the velocity fluctuations. Following [START_REF] Sipp | Characterization of noise amplifiers with global singular modes: the case of the leading-edge flat-plate boundary layer[END_REF], the singular values σ i (ω) are computed by generating the relevant operators in FreeFem++ and reformulating (3.17) into the following eigenvalue problem:

H * (ω)H(ω) φi (ω) = σ 2 i (ω) φi (ω). (3.18)
The largest eigenvalues of the Hermitian operator H * (ω)H(ω) are computed using the ARPACK library [START_REF] Lehoucq | Deflation techniques for an implicitly restarted Arnoldi iteration[END_REF] and the MUMPS parallel solver [START_REF] Amestoy | A fully asynchronous multifrontal solver using distributed dynamic shedding[END_REF]). Alternatively, one could use matrix-free methods by time-stepping the harmonically forced linearised Navier-Stokes equations (e.g. [START_REF] Monokrousos | Global threedimensional optimal disturbances in the blasius boundary-layer flow using time-steppers[END_REF]Gómez et al. 2016b). Since finite elements are used in this study, the Jacobian matrices are sparse and can be computed explicitly.

Reduced-order modelling of the fluctuation fields

In this section, we discuss how an efficient basis for the fluctuations can be educed from resolvent analysis by approximating the linear operator H(ω), if it is low-rank, or the nonlinear forcing f (ω). We also explain how time-resolved probe points are used to calibrate the amplitude and phase of the resolvent modes.

Approximation of the (linear) resolvent operator

The resolvent operator's singular values are said to be well separated if σ 1 (ω) σ 2 (ω) (similar to SPOD, the resolvent modes do not come in pairs for this flow), in which case the following approximation holds

H(ω) ≈ ψ1 (ω)σ 1 (ω) φ1 (ω).
(3.19)

(3.19) can be interpreted as an approximation of the (linear) resolvent operator and signifies that the first response mode ψ1 (ω) is sufficient to represent the velocity fluctuations at a frequency ω. In other words, it is possible to write

û(ω) ≈ ψ1 (ω)σ 1 (ω) φ * 1 (ω) f (ω) = ψ1 (ω)χ 1 (ω), (3.20)
where χ 1 is the complex weight of the resolvent mode obtained by projecting f onto the first basis vector for the nonlinear forcing φ1 and scaling by σ 1 . The projection φ * p (ω)f (ω), where p is an arbitrary forcing mode, is a measure of how effective the nonlinearity is in driving the velocity fluctuations by quantifying the angle between them. Therefore, the absolute magnitude of σ 1 alone does not determine what is observed in the flow although large values tend to coincide with energetic frequencies as will be seen in §7. The procedure for computing the complex weight χ 1 (ω) is described in §3.4.3. It should be noted that, in general, the Fourier mode û(ω) is not well defined unless the flow is truly periodic like cylinder flow. Hence, we refer to û(ω) and f (ω) as an arbitrary sample related to a specific bin.

Approximation of the (nonlinear) forcing

When there is an energetic frequency in the flow where the singular values are not well separated, the resolvent modes are no longer an efficient basis to represent the velocity fluctuations. An alternative approach is to approximate the nonlinear forcing f (ω) instead of the linear operator itself [START_REF] Rosenberg | The role of parasitic modes in nonlinear closure via the resolvent feedback loop[END_REF]

. (3.20) is modified into the following form û(ω) ≈ H(ω) f a (ω), (3.21)
where f a (ω) is the approximated nonlinear forcing that yields the parasitic mode. As stated in [START_REF] Mckeon | Experimental manipulation of wall turbulence: A systems approach[END_REF], the full expression for f (ω) is a convolution over all triadically consistent Fourier modes

f (ω c ) = ωa+ω b =ωc ωa,ω b =0 -∇ • (û(ω a )û T (ω b )). (3.22)
In this article, we will explain the conditions under which the sum in (3.22) can be severely truncated to a few key triadic interactions and for which types of flows this procedure is applicable.

Weighting the modes

Once a resolvent or parasitic mode for each frequency of interest has been computed using one of the two methods described above, its complex weight χ 1 (ω) is calibrated using a time-resolved probe ǔ(t, x 0 ) in a region of the flow where ψ1 (ω) is energetic as demonstrated by Gómez et al. (2016a) and [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF]. The signal is Fouriertransformed in time so that χ 1 (ω) can be computed in the frequency domain through the following expression

χ 1 (ω) = ǔ(ω, x 0 )/ ψ1 (ω, x 0 ).
(3.23) It should be noted that while ψ1 has been used in (3.23), the equation is still applicable to parasitic modes once their shapes have been computed using (3.21). After the complex amplitude of each mode has been determined, the linear superposition of all modes can be inverse Fourier-transformed to obtain a reduced-order model of the velocity fluctuations in the time domain, i.e.

u (t) = Nω j=1
ψ1 (ω j )χ 1 (ω j )e i(ωj t) .

(3.24)

3.5. Summary The resolvent-based modelling procedure is summarised in figure 2, which contains the relevant steps starting with the acquisition of the experimental measurements (top row). The mean flow is data-assimilated onto a FreeFem mesh with the correct boundary conditions and the time-resolved probe measurements are Fourier-transformed in time (row 2). The structure associated with frequency ω is obtained from resolvent analysis of the data-assimilated mean flow in one of two ways: (1) approximating the operator via the singular value decomposition if σ 1 σ 2 or (2) approximating the nonlinear forcing if the frequency is energetic but there is no separation of singular values. These two methods are demarcated by solid and dotted lines, respectively, in row 3 of figure 2. The resolvent modes are weighted from the Fourier-transformed probe points and inverse Fourier-transformed to obtain a model of the velocity fluctuations in time (bottom row).

Assimilated flow fields

In this section, the mean flows obtained from the time-averaged PIV are dataassimilated. Since data-assimilation attempts to reconstruct the optimal forcing f in the RANS equations, the spanwise component of ∇ × f , in addition to the assimilated mean fields, are compared to their experimental counterparts.

A0 Case

As mentioned in §2, the data below the airfoil are obstructed by the airfoil's shadow. To overcome this problem, the data above the centreline are reflected across, keeping the streamwise velocity constant and reversing the sign of the transverse velocity. Both the streamwise and transverse component of the data-assimilated fields are compared to their experimental counterparts in figure 3. The fields are nearly indistinguishable although the data-assimilated contours are slightly more smooth than those in the experiment.

The spanwise component of ∇×f , rather than f itself, is shown in figure 4 to eliminate the irrotational component which cannot be captured. The agreement between the experimental field, obtained after computing two gradients of the velocity fluctuations, and the data-assimilated field is good despite some minor discrepancies. One notable aspect about figure 4 is that ∇ × f for this flow closely resembles that of cylinder flow (see figures 4 and 5.9 of [START_REF] Foures | A data-assimilation method for Reynolds-averaged Navier-Stokes-driven mean flow reconstruction[END_REF] and Symon et al. (2018), respectively).

A10 Case

Data-assimilation of the A10 case is slightly more difficult since the flow is not symmetric with respect to the centreline. Two separate experiments are conducted at both α = -10 • and α = 10 • . In the first experiment, the data below the airfoil (suction side), are obstructed by the airfoil shadow so only the pressure side can be measured. In the second experiment, the pressure side is obstructed by the shadow, so only the suction side is measured. The final mean velocity field can be obtained by stitching together the two means such that the data around the entire airfoil are known. The stitching process, if imperfect, may result in non-physical measurements although the data-assimilation algorithm is able to compensate for this.

Data-assimilation is performed using the stitched mean velocity field, i.e. from measurements at α = -10 • and α = 10 • , and the resulting fields are compared to their experimental counterparts in figure 5. In addition to comparing the two components of the mean velocity field, the mean spanwise vorticity ω z is plotted in the bottom row of figure 5. Streamlines of the mean velocity field indicate the size and location of the recirculation bubble. The streamwise and transverse velocities match very well, particularly with respect to the reverse flow region. The separated flow region is accurately reproduced in the data-assimilated fields.

There are two regions of the flow field which are visibly improved by data-assimilation. The first is near the leading edge of the airfoil, where the flow separates and vortices are formed due to the rollup of the shear layer. This improvement is particularly visible in the vorticity plots since the negative contours of the data-assimilated field are continuous unlike in the experiment where there is a gap at x ≈ 0. in the outward movement of tracer particles in a direction perpendicular to the curved streamlines. This loss of seeding compromises measurements in the region close to the solid body.

The second area of improvement is at the trailing edge where there are imperfections associated with the stitching together of the two means. Contours of negative mean streamwise velocity protrude below the trailing edge which corresponds to non-physical behaviour. Data-assimilation is able to correct this problem and eliminate inaccuracies associated with the stitching process. These changes are more apparent by examining the difference between the assimilated and experimental mean velocities, i.e. ∆u = uu exp , which is illustrated in figure 6 for both velocity components. The two regions of disagreement, notably, coincide with the areas of improvement discussed earlier.

It should be emphasised that the benefit of data-assimilating fields is not only filling in missing data, but also in obtaining gradient quantities with less noise. This is particularly noticeable when comparing the spanwise component of ∇ × f in figure 7. Even though the experimental data do not have much noise to begin with, the data-assimilated field is more smooth overall.

Data-assimilation leads to improvement for both the A0 and A10 cases although the changes for A0 are fairly unremarkable. For the A10 case, on the other hand, there are features missing from the experimental mean that data-assimilation fills in, particularly at the shear layer originating near the leading edge. As will be seen in §6, this has major implications for the identification of linear amplification mechanisms in the flow.

SPOD results and the difference between an oscillator and an amplifier

In this section, SPOD is applied to the velocity fluctuations obtained from the raw PIV data. We also apply SPOD to the nonlinear term u • ∇u to investigate the most energetic structures associated with the nonlinear forcing field. [START_REF] Towne | Stochastic and nonlinear forcing of wavepackets in a Mach 0[END_REF] and . The SPOD eigenspectra of the A0 case (left) as a function of ω where red patches denote frequencies where there is significant separation between the first two eigenvalues. For the A10 case (right), the red patch denotes low-rank behaviour corresponding to the wake dynamics and the blue patch denotes low-rank behaviour corresponding to the shear layer dynamics.

Towne ( 2016) applied SPOD to the nonlinear forcing in a turbulent jet and found that the educed structures were primarily incoherent turbulent fluctuations and not nonlinear interactions between coherent structures. We will compare the results of SPOD applied to the velocity fluctuations and nonlinear forcing for both airfoils and highlight fundamental differences between oscillator and amplifier flows. These will have consequences on when we expect SPOD and resolvent analysis to coincide in §6 as well as our ability to model the flows in §7.

SPOD of the velocity fluctuations

The SPOD eigenspectra for A0 are plotted in figure 8 as a function of ω. The largest peak in λ 1 occurs at ω = 12.0, which corresponds to the shedding frequency of the airfoil. Red patches highlight frequencies at which there is significant separation between the first and second eigenvalues, the former of which exhibits sharp peaks at the shedding frequency and its harmonics.

The v-component of the first two SPOD modes is shown in figure 9 for the first three harmonics. In all cases, the first eigenvalue is at least one order of magnitude larger than the second eigenvalue. The first SPOD mode is significantly more structured than the second mode for all three harmonics in figure 9, an observation which is consistent with the fact that nearly all the energy at these frequencies is concentrated in the first mode. For frequencies where λ 1 ≈ λ 2 (not shown in the interest of brevity), the mode shapes are less structured for all eigenvalues.

Let us now consider the A10 case. The PIV data are subdivided into the same number of blocks containing the same number of snapshots as in the A0 case to allow for direct comparisons between both angles of attack. The SPOD eigenspectra are plotted in figure 8. It is immediately apparent that the SPOD eigenvalues are far less peaked than they were for the A0 case. Instead of observing a fundamental frequency and its harmonics, there are bands of frequencies where low-rank behaviour can be observed, a feature reminiscent of an amplifier flow. The first region, denoted by the red patch in figure 8, is situated around 4 < ω < 8 while the second one resides roughly between 10.5 < ω < 14.5.

The modes for A10 are plotted in figure 10 for the highest eigenvalues in the lowand high-frequency bands. Similar to the A0 case, the leading SPOD mode is highly structured and contains about one order of magnitude more energy than the second mode for ω = 6.01 and ω = 7.85. Most of the energy in these modes, furthermore, resides in the wake region, or x > 1, although there is some structure in the shear layer. The second SPOD mode, on the other hand, is more energetic in the shear layer although it is far less structured. The final frequency considered is ω = 12.47, which is roughly in the middle of the second band of frequencies identified earlier. The first mode is primarily in the shear layer while the second mode resembles a mix of both the shear layer and wake dynamics. Since SPOD takes into consideration the kinetic energy of the entire domain, the SPOD eigenvalues for frequencies corresponding to the shear layer assume lower values. They are dynamically important, nonetheless, since there is non-negligible separation between the first and second eigenvalue. They are also essential to model the fluctuations above the airfoil. The red patches in the A0 delineate the shedding frequency and its harmonics where there is separation between the first and second eigenvalues.

SPOD of the nonlinear forcing

We now utilise SPOD to characterise the structure of the nonlinear forcing, which is computed from the PIV snapshots since the data are sufficiently well-resolved in space and time. The SPOD for both (A0 & A10) are plotted together in figure 11 to underscore the contrast between the two cases. The eigenspectra for A0 are still peaked in the sense that the energy is high at the vortex shedding frequency and its harmonics. There is also a considerable gap between the first and second eigenvalues at these frequencies, similar to the trends observed for the velocity fluctuations. The spectra for the A10 case, on the other hand, are relatively flat and have no discernible peaks. The behaviour of the nonlinear forcing bears no resemblance to those of the velocity fluctuations as seen in figure 10.

The differences are even more evident when considering the mode shapes. For clarity, we will plot modes at the same frequencies as those in figures 9 and 10 since they are dynamically important. The modes for A0 are plotted in figure 12. It is clear that the nonlinear forcing for the first two harmonics is very structured, as the contours are almost perfectly symmetric (or anti-symmetric in the case of the second harmonic) with respect to the centreline. The structure corresponding to the third harmonic looks slightly more disordered, since it is not particularly energetic, but it can be remarked that the mode shape is approximately symmetric with respect to the centreline and has a well-defined wavelength which is approximately constant, i.e. not changing as a function of x.

The mode shapes for the A10 case are plotted in figure 12. The results are reminiscent of those from the jet flow study of [START_REF] Towne | Stochastic and nonlinear forcing of wavepackets in a Mach 0[END_REF] in that there is a lack of coherent structure. While there may be some notion of a wavelength associated with these modes, it would be far more difficult to quantify in comparison to the results in figure 12. There is also no noticeable difference between the spatial support of the lower frequency modes figure 13 and that of the higher frequency mode even though the choice of frequency has an impact on the spatial location and distribution of the velocity fluctuations. 5.3. Discussion and classification of the two cases SPOD of the velocity fluctuations reveals low-rank behaviour for both the A0 and A10 cases although there are some major differences. The SPOD eigenspectra of A0 are highly peaked at the vortex shedding frequency and its harmonics. The A0 case, thus, can be considered sparse in that there are few energetic frequencies in the flow. The SPOD eigenspectra for A10 are more broadband, with separation between the first and second eigenvalues occurring at two distinct ranges of frequencies making the flow less sparse, i.e. more frequencies are required to model the fluctuating flow field. Within the lower range, the SPOD modes are concentrated in the wake whereas within the higher range, they are concentrated in the shear layer. Based on these trends, the A0 case resembles an oscillator with fluctuations that are characterised by a well-defined frequency. The A10 case behaves more like an amplifier due to the shear layer instability although it can also be characterised as a weak oscillator due to its identifiable shedding frequency. Application of SPOD to the nonlinear forcing confirms this appraisal. In the case of A0, there is low-rank behaviour at all harmonics and the mode shapes are highly structured as they contain a well-defined wavelength and symmetry. This suggests they can be represented by a limited number of nonlinear interactions whereas in the A10 case, there is no apparent low-rank behaviour and the forcing field contains no coherent structure. These results are consistent with those of [START_REF] Rosenberg | The role of parasitic modes in nonlinear closure via the resolvent feedback loop[END_REF] and [START_REF] Towne | Stochastic and nonlinear forcing of wavepackets in a Mach 0[END_REF].

It is also in tune with the idea that the input to an amplifier is less important, as the linear dynamics have a preferential bias towards one particular structure.

In theory, this makes modelling of amplifier flows easier as the linear dynamics take precedence over nonlinearity, which needs to be taken into account for oscillator flows. Consequently, amplifier flows are more amenable to stochastic forcing of an approximated linear operator. Alternatively, the modelling of oscillator flows may benefit from approximating the nonlinear forcing due to the emergence of parasitic modes. These ideas are tested in §7.

Resolvent analysis of the data-assimilated means

For each angle of attack, two mean velocity fields are used as an input to resolvent analysis. The first is from the raw experimental mean, which is interpolated onto the same FreeFem mesh used in the data-assimilation algorithm. Due to the sparsity of data near the surface of the airfoil, the interpolated mean is not likely to satisfy the no-slip boundary condition, so this is manually enforced by setting the velocity at all points along the airfoil surface to zero. The points outside the PIV domain, where velocity vectors are not available, are set to uniform flow, i.e. u = (1, 0). This may seem like an oversimplification, but it has no discernible impact on the shape of the resolvent modes in the region where the experimental mean is known. The second mean is the dataassimilated mean, whose results are compared to those of the interpolated experimental mean in order to reinforce the benefits of data-assimilation before performing resolvent analysis. This section is divided into three parts: the results for the A0 and A10 cases are discussed in §6.1 and §6.2, respectively. In §6.3, the results from resolvent analysis are compared to those from SPOD and we remark on the different behaviours observed for an oscillator versus an amplifier flow.

A0 case: a single linear mechanism

The first three singular values of the resolvent operator are plotted in figure 14 using both the interpolated experimental and data-assimilated means. The trends for both cases are quite similar to cylinder flow (Symon et al. 2018) in that there is a distinct range of frequencies where the first singular value is an order of magnitude higher than the second singular value. Data-assimilating the mean velocity field has a negligible impact on the relative magnitudes of the singular values with respect to one another but it does increase the absolute value of the first singular value near the peak frequency. However, there is no impact on the mode shapes (see Symon 2018).

A10 case: multiple linear mechanisms

Unlike the A0 case, the singular values, illustrated in figure 15, of the resolvent operator for the A10 case depend strongly on the mean velocity field. The interpolated mean has a single peak at ω = 5.9 whereas the data-assimilated mean has two peaks at ω = 7.3 and ω = 15.4. The frequency range over which the resolvent is low-rank is also much wider for the data-assimilated mean.

The identification of a secondary peak can be attributed to filling in data in the shear layer near the leading edge. As seen in figure 16, which contains the most amplified resolvent mode at two frequencies corresponding to the peaks in figure 15 as well as one intermediate frequency (ω = 12.0), two separate linear mechanisms are identified. The first corresponds to shedding-or wake-type modes which occur at lower temporal frequencies (top two rows of figure 16). The second corresponds to a Kelvin-Helmholtz instability mechanism in the shear layer which occurs at higher temporal frequencies (bottom row of figure 16.) The absence of PIV vectors near the leading edge results in resolvent analysis failing to capture the second linear mechanism when the interpolated experimental mean is analysed instead of the data-assimilated mean. The mode shapes for intermediate frequencies have a mix of both linear mechanisms as seen for ω = 12 in figure 16. These results are consistent with [START_REF] Thomareis | Resolvent analysis of separated and attached flows around an airfoil at transitional Reynolds number[END_REF] and [START_REF] Yeh | Resolvent-analysis-based design of airfoil separation control[END_REF], who computed resolvent modes for a NACA 0012 airfoil at similar Reynolds numbers and angles of attack from time-averaged DNS and large-eddy simulation (LES) data, respectively. The latter study also identified two branches in the eigenspectrum of the LNS operator which were grouped into wake or shear layer modes.

6.3.

Comparison to SPOD Even though data-assimilation adds no new insight for the A0 case, there are several noteworthy differences between the results of resolvent analysis and SPOD for this flow. The most significant one is that resolvent analysis identifies a single peak at the shedding frequency and low-rank activity is limited to nearby frequencies only. There is no lowrank behaviour at higher harmonics whereas in SPOD, there is an order of magnitude separation between the first and second eigenvalue of the first three harmonics.

These observations can be explained by accounting for the nature of the nonlinear forcing. As in [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF], let us rewrite the cross-spectral density at an arbitrary frequency ω in terms of an expectation operator E(•):

Ŝ û û(ω) = E{û(x, ω)û * (x, ω)}.
(6.1) Noting that û(ω) = H(ω) f (ω), (6.1) can be rewritten as

Ŝ û û(ω) = E{H(ω) f (x, ω) f * (x, ω)H * (ω)}. (6.2)
Instead of expressing H(ω) in terms of its singular value decomposition, the resolvent operators are moved outside the expectation operator since they are deterministic to yield the following equation (6.3) where Ŝ f f (ω) = E{ f (x, ω) f * (x, ω)} is the cross-spectral density of the nonlinear forcing as in [START_REF] Towne | Resolvent-based estimation of space-time flow statistics[END_REF]. Ŝ û û(ω) may be low-rank in two scenarios. The first is when the resolvent operator is low-rank, which occurs when there is a linear amplification mechanism. This is observed at the shedding frequency of the A0 case, frequencies corresponding to wake modes in the A10 case, and frequencies of the shear layer modes in the A10 case as summarised in rows 1, 4 and 5 in table 1. Both SPOD and resolvent analysis identify low-rank behaviour, similar to the case of the turbulent jet in [START_REF] Schmidt | Spectral analysis of jet turbulence[END_REF] where the Kelvin-Helmholtz instability accounts for the low-rank nature of the cross-spectral density tensor.

Ŝ û û(ω) = H(ω) Ŝ f f (ω)H * (ω),
The second source of low-rank behaviour can originate from the cross-spectral density of the nonlinear forcing, or Ŝ f f (ω), being low-rank. Even if is nearly full-rank, as it is for higher harmonics of the shedding frequency, the left-hand side of (6.3) will be low-rank as long as at least one of the terms on the right-hand side is low-rank, i.e. rank( Ŝ û û(ω)) min rank(H(ω)), rank( Ŝ f f (ω)) .

(6.4)

For the A0 case, the nonlinear forcing is very structured at the shedding frequency and its harmonics; consequently, SPOD of the velocity fluctuations reveals large spectral gaps at these frequencies as seen in rows 1 and 2 of table 1. We theorise that the nonlinear forcing can be estimated by limited triadic interactions of resolvent modes to educe the structure at these higher harmonics, a hypothesis which is tested in the next section.

The SPOD eigenspectra of the nonlinear forcings in figures 12 and 13 verify that the cross-spectral density is low-rank at the shedding frequency and its harmonics in the A0 case and at no frequencies for the A10 The trends the notion that A0 behaves like an oscillator while A10 behaves more like an amplifier. It also suggests that resolvent modes and SPOD modes will deviate substantially for the higher harmonics in the A0 case since the resolvent is not low-rank and the forcing is anything but white noise in space and time. This will be seen in §7.1 for the second harmonic of the shedding frequency. The modes for the A10 case, on the other hand, can be expected to match at all frequencies where the the resolvent is low-rank since these coincide with the low-rank frequencies identified by SPOD.

We compare the dominant SPOD mode to the optimal resolvent mode in figure 17 to reinforce the notion that the first resolvent response mode is representative of the fluctuations in the flow when the resolvent is low-rank. One difference, which is particularly notable for the A0 shedding mode, is that the amplitude of the resolvent mode increases as a function of x. By comparing the final red contour in the top row of figure 17, we note that the amplitude of the SPOD mode is quite weak at x = 2.8, whereas it is very strong for the resolvent mode. A similar phenomenon can be observed for the A10 wake mode, which may result from A10 being a weak oscillator. The shear layer mode, on the other hand, is a very close match. Despite these discrepancies, the wavelengths of the modes are in complete agreement, signifying that the resolvent modes correctly predict the convection velocities of the structures in the flow.

Resolvent-based modelling of the velocity fluctuations

We now move on to determine how much of the full fluctuating field can be deduced starting from only a limited mean velocity field and data at a few probe points. We search for an efficient basis for a reduced-order model of the energetic coherent structures using the data-assimilated mean and very limited time-resolved probes (i.e. knowledge of the velocity at one or two points). This section is divided into three parts: in §7.1 we attempt to approximate the structure of the nonlinear forcing using limited triadic interactions at frequencies where the linear dynamics are not low-rank for the A0 case. The reducedorder models with and without these parasitic modes are compared. In §7.2, we consider the A10 case and examine whether just one or multiple probes are needed to account for the presence of two linear mechanisms in the flow. Finally, we discuss our results in the context of experimental limitations encountered in §7.3.

7.1. Parasitic modes and modelling of the A0 case

Approximation of the nonlinear forcing

To improve the basis for frequencies outside the range where the resolvent operator is low-rank, the nonlinear forcing is approximated using very limited triadic interactions. Let us consider two examples. The simplest is the second harmonic of the shedding frequency ω s , whose nonlinear forcing is approximated in exactly the same manner as that outlined in [START_REF] Rosenberg | The role of parasitic modes in nonlinear closure via the resolvent feedback loop[END_REF]

, i.e. f a (2ω s ) ∼ ψ1 (ω s ) • ∇ ψ1 (ω s ). (7.1)
Many other modes may contribute to the right-hand side of (7.1) although this approximation is adequate since the singular values associated with ψ1 (ω s ) are very large. Very low frequencies can also be recovered by postulating that they result from the beating of two highly amplified frequencies. The nonlinear forcing for the frequency ω = 3.70, for example, is

f (ω = 3.70) = - ωa+ω b =3.70 ωa,ω b =0 χ 1 (ω a )χ 1 (ω b )∇ • ( ψ1 (ω a ) ψT 1 (ω b )). (7.2)
Since the singular value associated with the shedding frequency at ω s = 12.0 is very large, it can be reasoned that the magnitude of χ 1 for this frequency will outweigh all other terms in (7.2). To obtain the correct frequency, therefore, it only needs to interact with the frequency ω = -8.30, where the resolvent operator is still approximately rank-1,

0 1 2 3 -0.2 0 0.2 0.4 0 1 2 3 -0.2 0 0.2 0.4 0 1 2 3 -0.2 0 0.2 0.4 0 1 2 3 -0.2 0 0.2 0.4 0 1 2 3 -0.2 0 0.2 0.4 0 1 2 3 -0.2 0 0.2 0.4
Figure 18. The optimal response mode (top row) compared to the parasitic mode (middle row) and the first SPOD mode (bottom row).

to produce a mode at the desired frequency ω = +3.70. (7.2), therefore, reduces to

f a (ω = 3.70) ∼ ψ1 (ω = 12.0)•∇ ψ1 (ω = -8.30)+ ψ1 (ω = -8.30)•∇ ψ1 (ω = 12.0), (7.3)
where the χ 1 have been dropped since χ 1 (ω = 12.0)χ 1 (ω = -8.3) appears before each term in (7.3) and, at this point, we are only interested in improving the basis for this frequency. If we were to include more triads in either (7.1) or (7.3), we would need to know the phase of the modes since χ 1 appears before each mode in (7.2). We argue here that the high amplifications associated with the shedding mode allow us to ignore the triadic interactions of suboptimal modes and frequencies which are less amplified.

The parasitic modes, or 'forced' resolvent modes (i.e. those obtained after rightmultiplying the resolvent operator by the approximated nonlinear forcing), are compared to the optimal response mode ψ1 from the singular value decomposition and the first SPOD mode in figure 18. Despite how simple the modelling for the nonlinear forcing is, the structure at all frequencies considered is a major improvement over the predicted structure from the first resolvent response mode. In particular, the parasitic mode obeys the correct symmetries while the optimal response mode resembles a stretched or compressed version of the vortex shedding mode.

Modelling preliminaries

Several candidate points to calibrate the complex weights of the resolvent modes are displayed in figure 19. They are referred to as points 1, 2, and 3, and are ordered by their streamwise position, i.e. point 1 is the furthest upstream. The power spectrum of P2 (P1 and P3 are very similar) for the transverse velocity component, plotted in figure 19, contains a large peak at = 12.0 and a secondary one at ω = 24.0. These are the first and second harmonics of the shedding frequency, respectively. There is another discernible peak at ω = 10.7 resulting from jitter in the shedding frequency. This is also present in the SPOD results for the leading eigenvalue, although it is less noticeable. The jitter does not influence the resolvent norm. The unweighted resolvent modes are first normalised using the standard energy inner product based on the two velocity components. Once the resolvent modes are weighted using the calibration point, they should reflect the true magnitude of the fluctuations. To make a fair comparison between the model and snapshots obtained from PIV, the experimental data are filtered to exclude temporal frequencies outside the range used in the model. The filtering involves Fourier-transforming the raw PIV data in time, truncating Fourier modes at unwanted frequencies, and inverse Fourier-transforming back into the time domain. Another possible way to compare the resolvent model to the experiment is to project the raw data onto the first SPOD mode at all frequencies of interest. Generally these would be where SPOD is approximately rank-1. Since the frequency resolution of the SPOD modes is only ∆ω ≈ 0.5, we choose to filter based on the first procedure.

Results and quantification of error

We begin by comparing the performance of each candidate point. The instantaneous error is computed between the model and PIV as follows

I v (t) = N i=1 v i (t) -v i P IV (t) 2 , (7.4)
where v denotes the transverse velocity component and N is the number of PIV vectors. I v (t) is integrated in time to obtain a global error metric

E v = 1 T T 0 I(t) 2 dt . (7.5)
E v is normalised by the total kinetic energy contained in the filtered PIV and the results for each point are summarised in table 2. Even though the performance of each point is good at some time instants, the overall performance for all three is relatively poor although P3 is the best. If two points are considered, the performance is nearly an order of magnitude better. Furthermore, points which are situated further downstream in the flow yield better results, which is a consequence of the resolvent modes failing to capture the correct streamwise decay of the fluctuations. The combination of points P2 and P3 yields the lowest amount of error, so the following discussion considers this choice only. The values computed in table 2 factor in all resolvent modes computed, i.e. both classical and parasitic. The first modelling attempt, however, only uses modes where the resolvent operator is nearly rank-1, i.e. 6.73 < ω < 15.7, henceforth referred to as When parasitic modes are included, the model is no longer constrained to be symmetric since the second harmonic and low frequency modes are antisymmetric with respect to the centreline. The flow is now modelled using 91 modes ranging from frequencies as low as ω = 2.72 to frequencies as high as ω = 39, in order to encompass the third harmonic. The resolvent-model is animated in time and compared to the filtered PIV in movies 1 and 2, respectively. Not all frequencies within this range are included since many of them are not energetic as indicated by the probe points. The model is compared to the filtered PIV in the bottom row of figure 20.

Because the modes outside the rank-1 range are not particularly energetic, the filtered PIV in the bottom row resembles the raw PIV snapshot in the top row fairly well (it contains 80% of the kinetic energy). The addition of parasitic modes does not dramatically improve the model, but it does capture more detailed features of the flow field. Further downstream, the negative velocity perturbations tilt upstream while the positive velocity perturbations tilt downstream and this is consistent with the PIV data. The discrepancy between the improved model and the filtered PIV is illustrated in figure 21. Small improvements can be seen in the region beyond x = 1.5. For this particular snapshot, the addition of parasitic modes reduces the error from 22% to 18%. 7.2. Two linear mechanisms and modelling of the A10 case

Selection of probe points

The challenge associated with the A0 case was to incorporate parasitic modes using nonlinear interactions. For the A10 case, the difficulty stems from being able to capture the two linear mechanisms present in the flow. Two candidate points are proposed in figure 22. P4 is located in the shear layer and P5 is in the airfoil wake. The power spectrum for P4 contains two broad peaks at frequencies similar to those identified by resolvent analysis and SPOD. This strengthens the argument that the resolvent tends to be lowrank at energetic frequencies as remarked by [START_REF] Moarref | Model-based scaling of the streamwise energy density in high-Reynolds-number turbulent channels[END_REF] for wall-bounded turbulent flows. The power spectrum for P5, on the other hand, contains one broad peak centred around the wake frequencies. It seems reasonable, therefore, to suggest that calibrating the modes using this point alone would not successfully model the fluctuations in the shear layer.

The performance using a single point or both points is summarised in table 3 using 66 modes in the range 3.14 < ω < 17.7. The reconstruction does very badly when only P4 is used as it overpredicts the fluctuations in the wake substantially. P5 is an improvement Table 3. Error for the A10 case as a function of point selection.

but it underpredicts the fluctuations in the shear layer. The model performs best when both points are used to calibrate the resolvent modes.

Modelling of the A10 velocity fluctuations

For the sake of brevity, only results for the use of both points are summarised in figure 23 as they help explain why a single point is not sufficient to obtain a good model of the flow. The raw PIV data are visualised in the top row. It is fairly evident that the flow is more complicated for this angle of attack as the presence of both linear mechanisms can be observed in relatively distinct regions of the flow. There is also an area immediately behind the trailing edge where they overlap leading to more complex flow structures.

Similar to the A0 case, the PIV data are frequency-filtered. We begin by considering 31 frequencies where roughly wake modes only are present, i.e. 3.15 < ω < 9.65, as seen in row 2 of figure 23. Even though there is some activity in the shear layer, it is small in comparison to the unfiltered case. The resolvent model is plotted alongside it. While the agreement between the wavelengths of the structures is excellent, the model overpredicts the intensity of the fluctuations as far upstream as x = 1.5.

The PIV data are next filtered to retain 31 frequencies corresponding to the shear Finally, the model is compared to the filtered PIV for 3.14 < ω < 17.7 (all 66 frequencies) in the bottom row of figure 23. Consistent with earlier results, the agreement is best in the shear layer while the performance of the model degrades further downstream. The model does capture some of the complex behaviour immediately behind the airfoil and the inclination of the fluctuations is correct overall. Despite poor performance further downstream, potential reasons for which are discussed in the next subsection, the resolvent model does quite well at reproducing the coherent structures of the flow.

Dynamics and convection velocities of the structures

Before discussing experimental limitations of the procedure, it is worth analysing the convection velocities of the structures in the resolvent model. To aid the reader, both the filtered PIV fluctuations (3.14 < ω < 17.7) and resolvent model are animated and uploaded as supplementary movies 3 and 4, respectively, to visualise the evolution of the flow in time. The movies contain 1000 snapshots and show that there is a coupling between the shear layer and wake modes in the region of the flow where the recirculation bubble is situated. Despite jitter in the flow's behaviour, the vortices which form in the shear layer meet those which emerge from the trailing edge and, in some cases, they merge. The convection velocity of the fluctuations, consequently, is relatively constant as a function of x although the perturbations begin to slow down further downstream. The reduced-order model from resolvent analysis is able to capture these key features despite poor predictions beyond x = 2.

Experimental limitations

Before we conclude, we discuss how the results could be improved and why this is difficult to achieve due to the limitations of the current experimental data sets.

The main problem associated with both cases is that the resolvent modes do not correctly model the streamwise decay of the fluctuations, a phenomenon also remarked upon by [START_REF] Rosenberg | The role of parasitic modes in nonlinear closure via the resolvent feedback loop[END_REF] for cylinder flow. In order to combat this problem, [START_REF] Rosenberg | The role of parasitic modes in nonlinear closure via the resolvent feedback loop[END_REF] approximated the nonlinear forcing for the shedding frequency using the triadic interaction between the shedding mode, predicted by the singular value decomposition, and the parasitic mode at the second harmonic, i.e. f a (ω s ) ∼ ψ1 (-ω s ) • ∇û a (2ω s ) + ûa (2ω s ) • ∇ ψ1 (-ω s ).

(7.6)

We attempted to do the same for the A0 and A10 cases since we have a good prediction of the fluctuation at the shedding frequency from the optimal resolvent mode and the second harmonic by approximating its nonlinear forcing. SPOD of the nonlinear forcing, furthermore, strongly suggests that it is necessary to account for this feedback since there is separation between the first two SPOD eigenvalues and the first mode is very structured. The approximated nonlinear forcing for the airfoil at the shedding frequency, however, did not improve the mode shape.

A possible reason for why this failed might be related to the enormous amplification at the shedding frequency due to the proximity of a stable eigenvalue to the neutral axis. Similar to cylinder flow, it is necessary to obtain a very accurate estimate of the nonlinear forcing in order to obtain the correct mode since the resolvent operator is so biased in favour of the first response mode. [START_REF] Lesshafft | Resolvent-based modeling of coherent wavepackets in a turbulent jet[END_REF] also remarked that spurious numerical noise can be sufficient to trigger the optimal mode and thus overwhelm suboptimal modes. The data-assimilated mean is not completely free of noise so the resolvent modes and their triadic interactions are not sufficiently accurate to alter the mode shape predicted by the resolvent at the shedding frequency. As such, an eddy viscosity model might improve the agreement between the SPOD and resolvent modes. Additionally, it would be worth investigating low Reynolds number DNS data, which are uncontaminated by noise, to understand the exact source of the deviations and how to alleviate these problems when using experimental data.

Since the nonlinear forcing is not particularly structured for the A10 case, it is also possible that three-dimensional effects may degrade resolvent-based predictions. Threedimensionality results from the formation of stall cells, the number of which is given by n s = AR/2.28 [START_REF] Weihs | Cellular patterns in poststall flow over unswept wings[END_REF], resulting in n s ≈ 2 for this setup. We observed earlier that the mean spanwise velocity component and its spanwise gradient are approximately zero. The mean flow, consequently, is nearly two-dimensional, hence the success of the data-assimilation (particularly with respect to the captured Reynolds stress gradients). In order to consider spanwise wavenumbers, however, we would need to obtain the spanwise-averaged mean velocity field (similar to [START_REF] Beneddine | Conditions for validity of mean flow stability analysis[END_REF], which is inaccessible from planar PIV. By considering spanwise wavenumbers other than k z = 0, which is more suitable for the region of the flow dominated by the Kelvin-Helmholtz instability, it may be possible to more accurately depict the far wake dynamics. The investigation of these proposed solutions is an area of current study.

Conclusions

The flows around an airfoil at two angles of attack, α = 0 • and α = 10 • , have been modelled from resolvent analysis using the procedure summarised in figure 2. The experimental measurements which were necessary included a crude mean velocity field obtained by PIV and no more than two probe measurements. SPOD of the velocity fluctuations indicated low-rank behaviour at the shedding frequency and its harmonics for the A0 case. The first eigenvalue of the cross-spectral density tensor was significantly higher than the second eigenvalue, resulting in large but narrow peaks. The A10 case, on the other hand, did not contain well-defined peaks although there was separation between the first two eigenvalues around two ranges of frequencies. The mode shapes within these ranges resembled wake and shear layer modes, respectively.

We also used SPOD to analyse the nonlinear forcing, which was directly computed from the PIV data. In the A0 case, the cross-spectral density of the nonlinear forcing exhibited similar trends to the velocity fluctuations in that a spectral gap existed at the shedding frequency and its harmonics. No such similarity existed for the A10 case as there was no low-rank behaviour at any frequency. The mode shapes observed for the A10 case, consequently, were highly unstructured, suggesting they did not emerge from a small number of modal interactions. As expected, the mode shapes for the A0 case were very structured for the shedding frequency and its harmonics, suggesting that the flow behaves like an oscillator with intrinsic dynamics which are insensitive to background noise. It followed that the A10 case behaved more like an amplifier.

In order to confirm this interpretation, the flows were analysed from a resolvent point of view. Before using the mean velocity field as an input, they were data-assimilated to decrease noise, fill in missing data, and correct artificial errors which resulted from stitching together multiple mean velocity fields (in the A10 case). The effect of dataassimilation was far more evident for the A10 case as it managed to identify two linear mechanisms which peaked at frequencies very similar to those identified by SPOD. When the mean velocity field was interpolated onto the CFD mesh, it was unable to identify the shear layer peak due to missing data near the separation point. Both the data-assimilated and interpolated means only identified a single linear mechanism for the A0 case at the shedding frequency, confirming our earlier interpretation that the A0 case behaved like an oscillator while the A10 case has characteristics of both types of flows.

By expanding the cross-spectral density tensor into the product of the resolvent operator, its Hermitian transpose, and the cross-spectral density tensor of the nonlinear forcing (similar to [START_REF] Towne | Spectral proper orthogonal decomposition and its relationship to dynamic mode decomposition and resolvent analysis[END_REF], we noted that low-rank behaviour may emerge when either the resolvent is low-rank (a linear amplification mechanism) or the crossspectral density tensor of the nonlinear forcing is low-rank, i.e. it is structured. The instances of low-rank behaviour due to either cause were summarised in table 1 for both airfoils. Notably in the A0 case, we educed the structure of the fluctuations at higher harmonics, where the resolvent operator is not low-rank, by considering limited triadic interactions of highly amplified resolvent modes.

When it came to resolvent-based modelling of the A0 flow, we no longer had to rely exclusively on the resolvent operator being low-rank to predict fluctuations. Indeed when the model was limited to rank-1 modes, we only recovered fluctuations which were symmetric with respect to the centreline since all rank-1 modes resembled the main instability mechanism. By adding parasitic modes, we were able to recover asymmetric fluctuations and capture finer details of the velocity fluctuations. Parasitic modes were not featured in the reconstruction of the A10 case. Instead, we determined that two probe points were necessary to capture both linear mechanisms in the flow and prevent overestimation of the fluctuations in the far wake region.

The combination of data-assimilation and resolvent analysis has shown its potential to fill in data from incomplete measurements and it is for this reason that they are a topic of ongoing study. The governing equations were exploited to obtain a mean velocity field which obeyed prescribed boundary conditions, removed experimental noise, and increased resolution. The performance of interpolation is clearly inferior as seen in the A10 case. Moreover, resolvent analysis produces a good basis which can be calibrated from limited time-resolved measurements.

It can also be concluded that the combination of SPOD and resolvent analysis can shed light on the nature of the nonlinear forcing in complex flows. In particular, they highlight when the nonlinear forcing is structured and the computation of parasitic modes in this study from experimental data suggests that this process can be leveraged for reducedorder modelling in the future. In cases where the nonlinear forcing is less structured, e.g. the A10 case, it could be possible to devise new strategies for predicting their statistics, e.g. [START_REF] Towne | Resolvent-based estimation of space-time flow statistics[END_REF], or modelling their effect on the large-scale structures, e.g. [START_REF] Illingworth | Estimating large-scale structures in wall turbulence using linear models[END_REF]. It may even suffice to rely on linear mechanisms only and treat the forcing as stochastic. The results from the paper also confirm the potential of resolvent analysis to predict the forces acting on a solid body, e.g. Gómez et al. (2016b). The A10 model, in particular, is quite promising in the vicinity of the airfoil so it would be possible to verify whether the predicted pressure field is valid by comparing the forces to those measured experimentally using a force sensor.

Figure 1 .

 1 Figure1. Sample vector plots of the instantaneous velocity for the A0 case (left) and A10 case (right). The mean flow has been subtracted and 1 out of every 2 vectors has been removed in x and y to facilitate visualisation of coherent structure.

Figure 2 .

 2 Figure 2. Steps to obtain the resolvent-based model of the velocity fluctuations starting from the experimental measurements, adapted from Symon (2018).

Figure 3 .

 3 Figure 3. Data-assimilation results for the A0 case. The assimilated velocity fields (right) are compared to the experiment (left). The streamwise velocity u is in the upper row while the transverse velocity v is in the lower row.

Figure 4 .

 4 Figure 4. The spanwise component of ∇ × f from the experiment (left) compared to the data-assimilation (right) in the A0 case.

Figure 5 .

 5 Figure 5. Data-assimilation results for the A10 case. The assimilated fields (right) are compared to the experiment (left): streamwise velocity u (top row), transverse velocity v (middle row) and spanwise vorticity ωz (bottom row).

  2. As noted by Raffel et al. (2018), there are difficulties associated with PIV for laminar flow around airfoils. The effects of strong centrifugal forces around the airfoil leading edge and strong shear result

Figure 6 .

 6 Figure 6. Residual discrepancy of the mean streamwise (left) and transverse (right) velocities between the experimental and data-assimilated mean velocity fields.

Figure 7 .

 7 Figure 7. The spanwise component of ∇ × f from the experiment (left) compared to the data-assimilation (right) in the A10 case.

  Figure8. The SPOD eigenspectra of the A0 case (left) as a function of ω where red patches denote frequencies where there is significant separation between the first two eigenvalues. For the A10 case (right), the red patch denotes low-rank behaviour corresponding to the wake dynamics and the blue patch denotes low-rank behaviour corresponding to the shear layer dynamics.

Figure 9 .Figure 10 .

 910 Figure 9. The first SPOD mode for the first three harmonics of the shedding frequency (left column) for the A0 case. The second SPOD mode for the first three harmonics of the shedding frequency (right column).

Figure 11 .

 11 Figure11. SPOD eigenspectra of the nonlinear forcing for the A0 (left) and A10 (right) cases. The red patches in the A0 delineate the shedding frequency and its harmonics where there is separation between the first and second eigenvalues.

Figure 12 .Figure 13 .

 1213 Figure 12. The first (left column) and second (right column) SPOD mode for the nonlinear forcing of the A0 case.

Figure 14 .

 14 Figure 14. The first three singular values of the resolvent operator plotted as a function of ω for the experimental (left) and data-assimilated (right) mean velocity fields of the A0 case.

Figure 15 .

 15 Figure 15. The first three singular values of the resolvent operator plotted as a function of ω for the experimental (left) and data-assimilated (right) mean velocity fields of the A10 case.

Figure 16 .

 16 Figure16. First optimal response mode ψ 1 for various frequencies using the A10 data-assimilated mean velocity field. The streamwise and transverse velocity components are plotted in the left and right columns, respectively.

Figure 17 .

 17 Figure17. A comparison between the dominant SPOD mode (left) and the optimal resolvent mode (right) at highly amplified frequencies.

Figure 19 .

 19 Figure19. Candidate probe points (left) and the power spectrum of v(t) for P2 (right).

Figure 20 .

 20 Figure 20. A0 results: (Top row) The unfiltered PIV snapshot. (Second row) Filtered PIV containing frequencies 6.73 < ω < 15.7 compared to the resolvent model. (Third row) Filtered PIV containing ω < 39 compared to the resolvent model.

Figure 21 .

 21 Figure 21. The discrepancy between the filtered PIV and the resolvent model using rank-1 frequencies only (left) and with the inclusion of parasitic modes (right).

Figure 22 .

 22 Figure 22. (Top) Probe points P4 and P5 in the shear layer and wake, respectively, and their power spectra (left and right) for the A10 case.

Figure 23 .

 23 Figure 23. A10 results: (Top row) The unfiltered PIV snapshot. (Second row) Filtered PIV containing frequencies 3.14 < ω < 9.88 compared to the resolvent model. (Third row) Filtered PIV containing 11.0 < ω < 17.7 compared to the resolvent model. (Fourth row) Filtered PIV containing 3.14 < ω < 17.7 compared to the resolvent model.

Table 1 .

 1 Properties of the cross-spectral density tensors and resolvent operators for the A0 and A10 cases.

	ω-range	Sûû(ω) (SPOD) H(ω) (Resolvent) S f f (ω)	modelling type
		A0 case (oscillator)
	1. near ωs	low-rank	low-rank	low-rank first resolvent mode
	2. near harmonics of ωs	low-rank	high-rank	low-rank	parasitic mode
	3. all other ω	high-rank	high-rank	high-rank	not modelled
		A10 case (amplifier)
	4. wake frequencies	low-rank	low-rank	high-rank first resolvent mode
	5. shear layer frequencies	low-rank	low-rank	high-rank first resolvent mode
	6. all other ω	high-rank	high-rank	high-rank	not modelled

Table 2 .

 2 Error for the A0 case as a function of point selection.

	1 Point Error 2 Points Error
	1	344% 1 + 2 165%
	2	988% 1 + 3	36%
	3	152% 2 + 3	32%

Acknowledgments

The support of the Army Research Office (ARO) through grant number W911NF-17-1-0306 and the Office of Naval Research (ONR) through grant number N00014-17-1-3022 are gratefully acknowledged.