Stéphane Perrard

Adrián Lozano-Durán

Marc Rabaud

Michael Benzaquen

Frédéric Moisy

Turbulent windprint on a liquid surface

Keywords: Free-surface flows, wall turbulence, wind waves

, the deformations that occur below the wave onset were named wrinkles. We provide here a theoretical framework for this regime, using the viscous response of a free liquid surface submitted to arbitrary normal and tangential interfacial stresses at its upper boundary. We relate the spatio-temporal spectrum of the surface deformations to that of the applied interfacial pressure and shear stress fluctuations. For that, we evaluate the spatio-temporal statistics of the turbulent forcing using Direct Numerical Simulation of a turbulent channel flow, assuming no coupling between the air and the liquid flows. Combining theory and numerical simulation, we obtain synthetic wrinkles fields that reproduce the experimental observations. We show that the wrinkles are a multi-scale superposition of random wakes generated by the turbulent fluctuations. They result mainly from the nearly isotropic pressure fluctuations generated in the boundary layer, rather than from the elongated shear stress fluctuations. The wrinkle regime described in this paper naturally arises as the viscous-saturated asymptotic of the inviscid growth theory of Phillips (1957). We finally discuss the possible relation between wrinkles and the onset of regular quasi-monochromatic waves at larger wind velocity. Experiments indicate that the onset of regular waves increases with liquid viscosity. Our theory suggests that regular waves are triggered when the wrinkle amplitude reaches a fraction of the viscous sublayer thickness. This implies that the turbulent fluctuations near the onset may play a key role in the triggering of exponential wave growth.

Introduction

A mirror like liquid surface is quite rare to observe in outdoor conditions. The smallest breeze already perturbs the surface of water well below the onset of wave formation, as first described by J.S. [START_REF] Russell | On waves[END_REF]. Since Russel's work, this weak deformation regime below the wave onset has been often reported [START_REF] Keulegan | Wind tides in small closed channels[END_REF][START_REF] Gottifredi | The growth of short waves on liquid surfaces under the action of a wind[END_REF][START_REF] Kahma | A laboratory study of the minimum wind speed for wind wave generation[END_REF][START_REF] Zhang | Capillary-gravity and capillary waves generated in a wind wave tank: Observations and theories[END_REF][START_REF] Caulliez | Drag of the water surface at very short fetches: Observations and modeling[END_REF]) but its precise spatiotemporal properties have been measured only recently by [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF] who named these deformations wrinkles. Characterised by randomly distributed, elongated structures in the streamwise direction, the apparent disordered aspect of the wrinkles were interpreted qualitatively as a signature of the air turbulence on the liquid surface.

An important feature of wrinkles is their dependence on the liquid viscosity ν . An empirical relation was inferred for the mean square amplitude [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF],

ζ 2 ∝ ν -1 u * 3 , (1.1)
where u * is the friction velocity of the air [START_REF] Schlichting | Boundary Layer Theory[END_REF] and ζ 2 is averaged over space and time. To the best of our knowledge, this empirical expression has not been explained theoretically.

Since wrinkles are systematically found even at moderate wind, they naturally represent a base state on which regular (quasi-monochromatic) waves propagating in the wind direction may grow. The influence of this initial perturbed surface state on the wave onset can be indirectly analysed by varying the liquid viscosity. Indeed, in the general understanding of the physical origin of the wind wave onset, the dependency in liquid viscosity is still lacking [START_REF] Sullivan | Dynamics of winds and currents coupled to surface waves[END_REF]. While it has been observed experimentally that the liquid viscosity influences the wave onset [START_REF] Francis | LXIX. Correspondence. Wave motions on a free oil surface[END_REF][START_REF] Kahma | A laboratory study of the minimum wind speed for wind wave generation[END_REF][START_REF] Veron | Experiments on the stability and transition of wind-driven water surfaces[END_REF][START_REF] Paquier | Viscosity effects in wind wave generation[END_REF], the explicit dependency has not been captured by models that include viscous effects [START_REF] Lindsay | The Kelvin-Helmholtz instability for a viscous interface[END_REF][START_REF] Funada | Viscous potential flow analysis of Kelvin-Helmholtz instability in a channel[END_REF][START_REF] Kim | Viscous effects on Kelvin-Helmholtz instability in a channel[END_REF].

In the literature of wind wave generation, two families of models can be identified. They involve either stability analysis of the mean wind profile or turbulent fluctuations. The branching goes back to the two seminal papers of [START_REF] Miles | On the generation of surface waves by shear flows[END_REF] and [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF] who proposed each a mechanism for wind wave generation. On the one hand, the stability analysis based on the mean wind profile in Miles models [START_REF] Miles | Surface-wave generation revisited[END_REF] or on the general Orr-Sommerfeld equation (see, e.g., Manneville 2010) ignores the turbulent fluctuations. Linear stability analysis predicts an onset above which the wave amplitude grows exponentially in time. A qualitative agreement with Miles theory has been obtained in laboratory experiments [START_REF] Kawai | Generation of initial wavelets by instability of a coupled shear flow and their evolution to wind waves[END_REF][START_REF] Veron | Experiments on the stability and transition of wind-driven water surfaces[END_REF]. However, a quantitative agreement of refined Miles models [START_REF] Janssen | The interaction of ocean waves and wind[END_REF]) with experimental data is still lacking both in laboratory conditions [START_REF] Plant | A relationship between wind stress and wave slope[END_REF][START_REF] Liberzon | Experimental study of the initial stages of wind waves' spatial evolution[END_REF] and in outdoor conditions [START_REF] Sullivan | Dynamics of winds and currents coupled to surface waves[END_REF]. On the other hand, Phillips, following [START_REF] Eckart | The generation of wind waves on a water surface[END_REF], analysed the effect of random pressure fluctuations at the surface of an inviscid liquid. He considered a resonance mechanism between the surface displacement ζ and the pressure fluctuations p. He derived an expression for the mean square surface displacement of the form

ζ 2 = 1 ρ 2 p 2 t 2 √ 2U c g , (1.2)
where p 2 is the mean square pressure fluctuation averaged over space, ρ the liquid density, U c the typical convection speed of the turbulent structures and g the acceleration of gravity. Phillips theory yields a linear growth for the wave energy. Contrary to Miles theory, Phillips theory has not been extensively tested. Experimentally, algebraic temporal growth for the surface deformation has been observed recently by [START_REF] Zavadasky | Water waves excited by near-impulsive wind forcing[END_REF]. Recent direct numerical simulations of two phase shear flows have also observed a regime of algebraic growth in time [START_REF] Lin | Direct numerical simulation of wind-wave generation processes[END_REF][START_REF] Zonta | Growth and spectra of gravity-capillary waves in countercurrent air/water turbulent flow[END_REF] that may be attributed to Phillips mechanism. However, a quantitative agreement with Phillips theory has not been reported to date. Phillips formalism may apply to the wrinkle regime, as the air flow is already turbulent for wind below the wave threshold, so that a resonance between pressure fluctuations and surface waves could already occur at low wind speed. Equation (1.2) provides a theoretical expression for ζ that depends on time and is independent of viscosity, whereas Eq. (1.1) provides an empirical relation for ζ that depends on viscosity and is independent of time. These two equations are therefore far apart elements of the wrinkle puzzle. The aim of this paper is to provide a theoretical framework for the wrinkle regime that reconciles Eqs. (1.1) and (1.2).

A quantitative description of the wrinkle regime, even though it does not involve any instability mechanism, is challenging for mainly two reasons. First, the formalism must describe the response of a viscous liquid to an arbitrary forcing both in time and in space. For an impulsive forcing of arbitrary shape one can follow the approach of [START_REF] Miles | The Cauchy-Poisson problem for a viscous liquid[END_REF], who revisited the Cauchy-Poisson problem [START_REF] Lamb | Hydrodynamics[END_REF] for a viscous liquid. The response to a continuous perturbation in time can also be computed using the same formalism by linear superposition. However, the Fourier-Laplace transform formalism used by Miles limits the analytical feasibility to asymptotic solutions. For the specific case of a pressure source travelling at constant velocity, the wave pattern generated at the surface of an inviscid liquid was computed by [START_REF] Havelock | Wave resistance: Some cases of three-dimensional fluid motion[END_REF]. This classical problem, which provides a simplified description of the far-field wake behind a ship, has been recently revisited, for an inviscid [START_REF] Raphaël | Capillary gravity waves caused by a moving disturbance: wave resistance[END_REF][START_REF] Rabaud | Ship wakes: Kelvin or Mach angle?[END_REF][START_REF] Darmon | Kelvin wake pattern at large Froude numbers[END_REF] or a viscous [START_REF] Richard | Capillary-gravity waves: The effect of viscosity on the wave resistance[END_REF]) fluid. To the best of our knowledge, no such Havelock-like formulation is yet available for an arbitrary forcing pattern both in time and in space. Moreover, the Havelock formulation applies only for a pressure disturbance, whereas both pressure and shear stress act on the surface of a viscous liquid.

The second main difficulty arises from the modelling of the turbulence in the air boundary layer. Of particular interest for the wave generation problem is the slow dynamics and the long-range correlations of the pressure and shear stress fluctuations in the boundary layer. Their main statistical properties (characteristic size, mean convection velocity) have been measured since the 60's, e.g. by [START_REF] Willmarth | Measurements of the fluctuating pressure at the wall beneath a thick turbulent boundary layer[END_REF] and [START_REF] Corcos | The structure of the turbulent pressure field in boundary-layer flows[END_REF]; see [START_REF] Robinson | Coherent motions in the turbulent boundary layer[END_REF] for a comprehensive review prior to the development of numerical simulations. These early quantitative measurements, obtained from singlepoint probes, however, could not provide a full spatio-temporal characterisation of the turbulent fluctuations in a boundary layer.

In recent years, in-depth knowledge has been gained from Direct Numerical Simulations (DNS), both for developing turbulent boundary layers and for fully developed turbulent channel flows [START_REF] Choi | On the space-time characteristics of wall-pressure fluctuations[END_REF][START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Re = 590[END_REF][START_REF] Jimenez | The large-scale dynamics of near-wall turbulence[END_REF][START_REF] Jimenez | Turbulent fluctuations above the buffer layer of wall-bounded flows[END_REF]. The maximum turbulent Reynolds number reached in the most advanced simulations, Re τ 4000 -8000 (Lozano-Durán & Jiménez 2014; [START_REF] Lee | Direct numerical simulation of a turbulent channel flow up to Reτ ≈ 5200[END_REF][START_REF] Yamamoto | Numerical evidence of logarithmic regions in channel flow at Reτ = 8000[END_REF], is comfortably larger than the relevant values for the windwave generation problem (Re τ 100-1000). Such data are highly valuable for the study of the wrinkle regime, as the air turbulence can be considered as essentially unaffected by the wave motion. Indeed, the typical surface displacement amplitude is much smaller than the viscous sublayer thickness in the air flow. This is precisely the regime explored in this paper. The situation is different for larger amplitude, close to the wave onset, for which a feedback of the wave motion on the air turbulence is expected, requiring a full two-phase flow approach. This approach, much more demanding computationally speaking, has been investigated recently [START_REF] Belcher | Turbulent flow over hills and waves[END_REF][START_REF] Kudryavtsev | Coupled dynamics of short waves and the airflow over long surface waves[END_REF][START_REF] Lin | Direct numerical simulation of wind-wave generation processes[END_REF][START_REF] Druzhinin | Direct numerical simulation of a turbulent wind over a wavy water surface[END_REF][START_REF] Kudryavtsev | Impact of wind waves on the air-sea fluxes: A coupled model[END_REF][START_REF] Sajjadi | Wave motion induced by turbulent shear flows over growing Stokes waves[END_REF]. Because of the high computational cost, the range of physical parameters covered by these studies remains limited, and no scaling relation has yet been obtained for the wave onset.

An important feature in the generation of surface waves by the wind is the unavoidable presence of a mean advecting current in the liquid, which may significantly affect the wave dynamics [START_REF] Peregrine | Interaction of water waves and currents[END_REF][START_REF] Kirby | Surface waves on vertically sheared flows: Approximate dispersion relations[END_REF][START_REF] Ellingsen | Approximate dispersion relations for waves on arbitrary shear flows[END_REF]. The growth and saturation of a drift velocity in the liquid and wave growth are two closely intertwined processes [START_REF] Banner | Tangential stress beneath wind-driven air-water interfaces[END_REF][START_REF] Melville | Laboratory measurements of the generation and evolution of langmuir circulations[END_REF][START_REF] Veron | Experiments on the stability and transition of wind-driven water surfaces[END_REF], and the stability of the flow depends on the resulting velocity profiles both in the air and in the liquid [START_REF] Miles | Surface-wave generation revisited[END_REF]. However, in the wrinkle regime explored by [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF], the drift velocity is bounded to a few cm/s by the finite channel depth and the large viscosity of the fluid. Such drift velocity remains small compared to both air velocity and phase velocity of the wrinkles, and is neglected in the present study.

In this paper we focus on a quantitative description of the wrinkle regime, below the wind wave onset. For that purpose, we combine analytical calculations of the viscous response of the liquid to an arbitrary forcing in a statistically stationary state with DNS data of a turbulent channel flow. Based on careful dimensional analysis and experimental data of the wrinkle regime, we show that below the wave threshold the evolution of the turbulent air boundary layer can be decoupled from the liquid response. We thus circumvent the difficulty to simulate a full turbulent two-phase flow by considering a turbulent channel flow with rigid walls and no-slip boundary conditions. Doing so, we greatly simplify the numerical setup to focus on the linear, passive response of the liquid phase. We show that the wrinkle statistics computed from our model are in good agreement with the experimental data of [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF][START_REF] Paquier | Viscosity effects in wind wave generation[END_REF]. Finally, we show that the onset of regular waves may correspond to the breakdown of the regime of linear passive response of the liquid surface, and propose an empirical criterion based on surface roughness originating from the wrinkle amplitude.

The paper is organized as follows. Section 2 provides a general dimensional analysis of the surface deformation generated by a turbulent boundary layer. In section 3 we derive the key equation of the paper, that establishes the link between the Fourier components of the normal and tangential shear stresses applied at the air-liquid interface and the Fourier components of the surface displacement. Section 4 presents the direct numerical simulations used to compute the normal and tangential stresses. Section 5 combines the equation for the surface displacement in Fourier space and the output of the DNS to compute the wrinkle properties. It provides a quantitative comparison with the experiments of [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF]. Section 6 finally bridges the gap between the wrinkle regime governed by viscosity and the inviscid resonant theory of [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF].

Wrinkle regime: dimensional analysis and experimental set up

Dimensional Analysis

We first discuss here the dependence of the surface displacement ζ on the relevant physical parameters using dimensional analysis. In the statistically steady state, the characteristic amplitude of the surface displacement averaged over space and time ζ rms = (ζ 2 ) 1/2 in response to a turbulent wind is expected to depend on numerous parameters that characterise both the turbulent flow in the air and the liquid properties. The geometry and the relevant parameters of both phases are sketched in figure 1. The turbulent air boundary layer is characterised by the air density ρ a , the kinematic viscosity ν a , the friction velocity u * = τ a /ρ a (where τ a is the mean shear stress at the interface), and the boundary layer thickness δ. In a developing boundary layer, the thickness δ is function of the streamwise distance x, usually called fetch in the wind-wave generation problem, whereas it is constant in a fully developed channel flow. The spatial variation of a developing boundary layer is usually small (dδ/dx 1), so we can simply consider the boundary layer thickness δ as constant. The liquid flow is characterised by the liquid

ζ(r, t) z x g Figure 1.
Sketch of the velocity profile of a turbulent wind blowing on a viscous liquid. The turbulent boundary layer is characterised by the outer layer thickness δ and the viscous sublayer thickness δν = νa/u * . The friction velocity u * is defined from the mean shear stress τa = ρau * 2 at the wall. The liquid response may depend on the air/liquid densities ρa, ρ , the kinematic viscosities νa, ν , as well as δ, u * , the gravity g and the surface tension γ. density ρ , kinematic viscosity ν , acceleration of gravity g, surface tension γ, and liquid depth h.

Under these hypotheses, the amplitude of the surface displacement writes ζ rms = f (ρ a , ρ , ν a , ν , δ, u * , g, γ, h). According to Buckingham's π theorem, the dimensionless surface displacement ζ rms /δ can be expressed as a function of six independent dimensionless numbers. Choosing a set of dimensionless numbers that decouple the influence of the friction velocity u * and of the length scale δ, it writes

ζ rms δ = f 1 ρ a ρ , gδ 3 ν 2 , u * 3 gν , Re δ , Bo, h δ , (2.1)
where f 1 is a dimensionless function, Re δ = u * δ/ν a is the turbulent Reynolds number characterising the boundary layer, and Bo = δ/ c is the Bond number (with c = γ/ρ g the capillary length). The relative depth h/δ becomes relevant for surface wave propagation in the shallow water regime. We are mostly interested here in the the deep water regime (one has h/δ 1.2 in the experiment, see section 2.2), so the importance of this parameter is marginal in the following. The function f 1 can be further specified using additional physical arguments. In the static case and without surface tension, the density ratio ρ a /ρ sets the surface displacement amplitude as can be seen from a simple pressure balance: the gravity pressure scales as ρ gζ and the pressure fluctuations in the air phase scales as ρ a u * 2 . Adding the surface tension introduces a dependency in Bond number but it does not modify the scaling of ζ rms in ρ a /ρ . In the dynamical case with negligible effect of gravity in the air phase, and in the limit of linear equation of motion, the displacement still scales as ρ a /ρ as can be seen from the continuity of normal stresses at liquid-air interface (a proper justification is given in Section 3). It implies

ζ rms δ = ρ a ρ f 2 gδ 3 ν 2 , u * 3 gν , Re δ , Bo, h δ . (2.2)
The dimensionless number gδ 3 /ν 2 characterises the surface response of a viscous liquid to an initial perturbation, as analysed by [START_REF] Miles | The Cauchy-Poisson problem for a viscous liquid[END_REF]. The dependency in gδ 3 /ν 2 can be neglected for the following reason. We first rewrite gδ 3 /ν 2 as (δ/ ν ) 3 , with ν = g -1/3 ν 2/3 the viscous length. This length ν was identified by [START_REF] Miles | The Cauchy-Poisson problem for a viscous liquid[END_REF] as the relevant length scale to classify the surface deformation regimes. For a gravity wave of wavenumber δ -1 and angular frequency ω = g/δ, the dimensionless damping factor defined as θ = ν /(ωδ 2 ) is given by θ = ( ν /δ) 3/2 . The separation between the propagating wave regime and the over-damped regime occurs at a finite value of θ, namely θ c = 1.31 [START_REF] Leblond | The viscous damping of capillary-gravity waves[END_REF]. Hence for gravity waves ν separates the regime of propagating waves (θ < θ c ) from the over-damped regime (θ > θ c ). In the range of liquid viscosity for which wrinkles are observed, ν 1 -10 3 ν water , the viscous length ν is in the range 50 µm -5 mm. Natural and laboratory flows characterised by a forcing scale δ >1 cm and for moderate viscosity ν < 10 3 ν water therefore fall in the regime δ ν in which gravity dominates over viscous diffusion. Surface deformations will therefore not be significantly diffused horizontally by viscous effects, and we may thus neglect the influence of gδ 3 /ν 2 on the surface deformation.

Even though the dimensionless number gδ 3 /ν 2 is non relevant in the propagating regime, viscosity cannot be neglected because its cumulative effect eventually balances the input forcing [START_REF] Miles | The Cauchy-Poisson problem for a viscous liquid[END_REF]). The displacement ζ rms may therefore depend on u * 3 /(gν ), and one is left with

ζ rms δ = ρ a ρ f 3 u * 3 gν , Re δ , Bo, h δ . (2.3)
Finally, from a balance between the energy flux from the turbulent boundary layer in the air and the viscous loss in the liquid, it is possible to infer the dependence of ζ in the liquid viscosity ν . Neglecting surface tension and finite depth effects (large Bo and large h/δ), the potential energy e per unit surface of a deformation of amplitude ζ is of order ρ gζ 2 . If this fluctuation is in the wave regime (weak viscous attenuation), its kinetic energy per unit surface is also of order e. We assume here that the liquid surface is mostly sensitive to the largest scales of the turbulent fluctuations, governed by the boundary layer thickness δ (we provide support to this key assumption in section 5). Consider a vertical velocity fluctuation of order u * over horizontal extent of order δ, corresponding to a pressure fluctuation of order ρ a u * 2 , pushing or sucking the liquid surface at a velocity ζ. Conservation of vertical momentum in this inelastic process implies ρ a u * ρ ζ. The power per unit surface transferred to the liquid, ρ a u * 2 ζ, therefore writes (ρ 2 a /ρ )u * 3 . In a statistically steady state, this power must be balanced by the energy loss by viscous diffusion in the liquid, e/τ , with τ δ 2 /ν the viscous time scale, yielding ζ δ(ρ a /ρ )(u * 3 /gν ) 1/2 . This suggests writing Eq. (2.3) in the form †

ζ rms δ = ρ a ρ u * 3 gν 1/2 f 4 Re δ , Bo, h δ .
(2.4)

The aim of this paper is to provide a solid mathematical ground to this scaling using viscous surface waves theory and an appropriate evaluation of the surface response in Fourier space following a route similar to [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF].

In the following, the full dependency of ζ/δ with respect to the Bond number and liquid depth is not explored, as the experimental data on the wrinkle regime are available only for one value of Bo and h/δ (see section 2.2 for details). The remaining dependency of ζ/δ in Re δ is subtle and rises the question of the relevant forcing scale in this problem. The above qualitative argument assumed that the forcing acts at the scale of the boundarylayer thickness δ, whereas the spectrum of the stress fluctuations in a turbulent boundary layer spreads in the range from δ down to the viscous sublayer thickness δ ν = Re -1 δ δ. The non-trivial dependence in Re δ therefore depends on the exact forcing spectrum and the nature of the surface response, and will be characterised in section 5 using DNS data.

Re δ = u * δ/νa 10 2 -10 3 Bo = δ/ c 14 h/δ 1.2 ρa/ρ 0.9×10 -3 -1.2×10 -3 u * 3 /gν 2 × 10 -2 -2 × 10 2 gδ 3 /ν 2 8 × 10 2 -3 × 10 8
Table 1. Set of non-dimensional numbers in the experiments [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF][START_REF] Paquier | Viscosity effects in wind wave generation[END_REF].

The key result of the paper is that, in the range of Re δ relevant for the experiments, the function f 4 is essentially constant, which is consistent with the empirical scaling (1.1).

Experimental details

Although the theory derived in this paper is general, quantitative comparison in the following is provided with the only available experimental data of [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF][START_REF] Paquier | Viscosity effects in wind wave generation[END_REF]. We briefly provide here some details about the experiment, and summarise the relevant non-dimensional numbers in table 1.

The experimental set up consists in a rectangular tank filled with liquid, fitted to the bottom of a wind-tunnel of rectangular cross-section. The tank is of length L = 1.5 m, width 296 mm, and depth h = 35 mm, and the channel height is 105 mm. Air (density ρ a 1.2 kg m -3 , kinematic viscosity ν a 15 × 10 -6 m 2 s -1 ) is injected at a mean velocity U a in the range 1 -10 m s -1 . The velocity profile in the air is close to that of a classical turbulent boundary layer developing over a no-slip flat wall, at least in the wrinkle regime. The boundary-layer thickness δ, defined as the distance from the surface at which the mean velocity is 0.99U a , is δ 13 mm at x = 0, and increases linearly along the tank, with dδ/dx 0.02. At the fetch x at which the measurements are performed, the local boundary layer thickness is δ 30 mm. The friction velocity u * , determined from measurement of the surface drift velocity using stress continuity, is approximately u * 0.05U a .

The liquid viscosity is varied in a wide range, ν = 0.9 -560 × 10 -6 m 2 s -1 , using mixtures of glycerol and water (for low ν ) or glucose syrup and water (for large ν ). The liquid density ρ of the mixtures is in the range 1.0 -1.36 × 10 3 kg m -3 . The surface tension is γ 60 mN m -1 , and the capillary length c = γ/ρ g is approximately 2.2 mm.

For a laminar flow in the liquid, the drift velocity at the surface is given by the continuity of the shear stress, namely U s = ρ a u * 2 h/(4ρ ν ). For liquid viscosity ν > 4ν water , the maximum drift velocity in the wrinkle regime is U s 10 cm/s (for a wind velocity U a 4.5 m/s), for which the flow in the liquid remains essentially laminar (U s h/ν < 10 3 ). For the largest liquid viscosity, ν 600ν water , the surface drift does not exceed 1 mm/s even at the largest wind velocity, and can be safely neglected. However, for water and liquids of viscosity up to 4ν water , the drift velocity is significant and the flow in the liquid is no longer laminar. For the sake of simplicity, the effects of the drift velocity and of turbulence in the liquid are not considered in the present paper.

The instantaneous surface deformation fields ζ(r, t) are measured using Free-Surface Synthetic Schlieren [START_REF] Moisy | A synthetic schlieren method for the measurement of the topography of a liquid interface[END_REF]. This optical method is based on the analysis of the refracted image of a pattern visualized through the interface. The field of view is 390 × 280 mm. The horizontal resolution is 3 mm, and the vertical resolution is 0.6 µm. The rms wrinkle amplitude is in the range 1-10 µm (ζ rms /δ 3 × 10 -5 -3 × 10 -4 ).

Derivation of the spectral theory

In this section we derive an expression in Fourier space relating the surface displacement to an arbitrary pressure and shear stress fields applied at a liquid interface from the upper gas phase.

Assumptions

The calculation is made under the following assumptions:

(a) The gas density is small compared to the liquid density, ρ a ρ . (b) The slope of the surface displacement stays small at all time (|∇ζ| 1). (c) The flow in the liquid is laminar. (d) The drift velocity in the liquid is negligible compared to the convection speed of the turbulent structures and the typical phase velocity of the wrinkles.

(e) The liquid layer is of infinite depth. Assumptions (a)-(c) are fundamental hypotheses of our theoretical approach. The drift current (assumption (d)) strongly affects the propagation of surface wave in oceanographic context [START_REF] Peregrine | Interaction of water waves and currents[END_REF], and correction to the dispersion relation has to be taken into account [START_REF] Kirby | Surface waves on vertically sheared flows: Approximate dispersion relations[END_REF][START_REF] Ellingsen | Approximate dispersion relations for waves on arbitrary shear flows[END_REF]. However, in our experimental range of parameters (see section 2.2), the drift current stays within few per cent of the convection speed of the turbulent structures and the typical phase velocity of the wrinkles, so it can be neglected. The finite depth correction (e) will be included in the limit of bulk-dominated dissipation in section 3.4. Two other assumptions shall also be used later but are not required for the following main derivation:

(f) The surface displacement falls into the propagative wave regime, i.e. the horizontal sizes of the surface displacement are larger than the viscous length ν = g -1/3 ν 2/3 .

(g) The surface displacement ζ is small compared to the viscous sublayer δ ν of the turbulent boundary layer in the air.

Assumption (f) is comfortably satisfied experimentally in the wrinkle regime for usual values of liquid viscosity, ν < 1000ν water . Assumption (g) will be useful in section 4 to model the pressure and shear applied at the liquid interface by the gas layer.

Formulation

We consider a liquid layer submitted to a normal stress field N (x, y, t) and a shear stress T (x, y, t) = T x e x + T y e y applied at its upper surface in z = ζ(x, y, t). The linearized Navier-Stokes equation for the velocity v = v x e x + v y e y + v z e z in the liquid reads

∂ t v = - 1 ρ ∇p + g + ν ∆v, (3.1)
where p is the pressure in the liquid phase. We have ∇ • v = 0 and the boundary condition for an infinite liquid depth reads

lim z→-∞ v = 0. (3.2)
At the interface in z = ζ, the continuity condition of the normal stress can be approximated by the pressure in z = 0 using p (x, y, z = ζ) = p 0 -ρ gζ, where p 0 = p (x, y, z = 0). The continuity condition for the normal stress in z = 0 yields

p 0 -ρ gζ -2ρ ν (∂ z v z ) z=0 + γ∆ (x,y) ζ = N, (3.3)
where ∆ (x,y) is the 2D Laplacian. The continuity condition of the tangential stress writes

ρ ν ∂ z v x + ∂ x v z z=0 = T x (3.4a) ρ ν ∂ z v y + ∂ y v z z=0 = T y .
(3.4b)

Statistically steady and homogeneous regime

We now look for statistically homogeneous and stationary solutions. We introduce the space-time Fourier transform F,

ζ(k, ω) = F{ζ(r, t)} = d 2 r dt ζ(r, t)e -i(k•r-ωt) (3.5a) ζ(r, t) = F -1 { ζ(k, ω)} = (2π) -3 d 2 k dω ζ(k, ω)e i(k•r-ωt) , (3.5b)
and similarly for N (k, ω) and T (k, ω), with r = xe x + ye y the horizontal position, k = k x e x + k y e y the horizontal wave vector and ω the angular frequency (k and ω are real). From Eq. (3.1) we obtain ∆p = 0 (3.6)

∂ t Ω -ν ∆Ω = 0, (3.7)
where we have introduced the vorticity Ω = ∇×v. The boundary conditions for Ω writes

Ω(r, z = 0, t) = Ω 0 (r, t) (3.8) lim z→-∞ Ω(r, z, t) = 0, (3.9)
where Ω 0 (r, t) is the vorticity at the surface. Eqs. (3.6) and (3.7) can be solved in Fourier space, (3.11) where

p (r, z, t) = (2π) -3 d 2 k dω p0 (k, ω)e i(k•r-ωt) e kz (3.10) Ω(r, z, t) = (2π) -3 d 2 k dω Ω0 (k, ω)e i(k•r-ωt) e mz ,
m 2 = k 2 -iω/ν , (3.12) with k = (k 2 x + k 2 y ) 1/2
and {p 0 , Ω0 } = F{p 0 , Ω 0 } being the Fourier transforms of the pressure and vorticity in z = 0. Eq. (3.10) shows that an applied pressure patch of typical wavenumber k and frequency ω penetrates the liquid layer over a depth k -1 , while Eq. (3.11) shows that the applied shear stress and hence vorticity penetrate over a depth |m| -1 . These two penetration depths are similar for very viscous fluids, whereas in the limit of low viscosity vorticity remains confined in a thin boundary layer of thickness

|m| -1 ν /ω. Rewriting Eq. (3.1) as ∂ t v = -∇p /ρ -ν ∇ × Ω,
v(r, z, t) = (2π) -3 d 2 k dω q ρ iω p0 e kz + ν iω κ × Ω0 e mz e i(k•r-ωt) , (3.13)
where q = (ik x , ik y , k) and κ = (ik x , ik y , m). The functions p0 and Ω0 are given by the stress boundary conditions in z = 0. Combining Eqs. (3.3), (3.4) and (3.13) yields

1 - 2k 2 ν iω k ρ p0 - 2ν mk iω ν Bz -g k ζ = k N ρ (3.14) ν Bz + 2ν F{(∂ zz v z ) z=0 } = - ik • T ρ , ( 3
ω 2 -g k + 4iν ωk 2 + 4ν 2 k 3 (m -k) ζ = k N ρ + m -k m + k ik • T ρ .
(3.17)

We obtain a damped wave equation forced by the normal and tangential stresses applied at the interface. It describes the response of a viscous liquid forced by arbitrary normal and tangential stress fields under the assumptions (a) to (e), in statistically steady and homogeneous configurations. Setting N =0 and T = 0 yields the usual dispersion relation for gravity-capillary waves with viscosity [START_REF] Lamb | Hydrodynamics[END_REF]. In a very viscous fluid (m k), the effect of the shear stress on the flow vanishes, whereas for a fluid of low viscosity (|m| k), the effects of pressure and shear stress are comparable. Note here the specificity of the limit ν → 0: although (mk)/(m + k) → 1, only pressure can generate waves since the inviscid limit implies T = 0 by tangential stress continuity.

Eq. (3.17) can be simplified using the small viscosity approximation (assumption (f)). The boundary layer thickness associated to the vertical diffusion of a shear stress patch applied during a time ω -1 is δ = ν /ω. For a gravity wave of frequency ω = √ gk we have δ = 3/4 ν k -1/4 , where ν = (ν 2 /g) 1/3 is the viscous length introduced in section 2.1. In practice, δ lies in the range 0.2 -7 mm for a typical wrinkle wavelength 2π/k ≈ 100 mm. The thin boundary layer approximation (f) is therefore fulfilled. In this limit, we have |m| k, and only the first order in ν contributes in Eq. (3.17), yielding

ζ(k, ω) = 1 ρ k N + ik • T ω 2 -g k + 4iν ωk 2 . ( 3
ζ(r, t) = 1 (2π) 3 1 ρ d 2 k dω k N + ik • T ω 2 -g k + 4iν ωk 2 . (3.19)
In the specific case of T (k, ω) = 0 and N (k, ω) in the form δ(ω -U c k x ) N (k), with U c the convection velocity of a rigid pressure source of Fourier transform N (k), this equation reduces to the classical Havelock integral used to describe the far-field wake of a ship [START_REF] Havelock | Wave resistance: Some cases of three-dimensional fluid motion[END_REF][START_REF] Raphaël | Capillary gravity waves caused by a moving disturbance: wave resistance[END_REF]. with the source term defined as The real part of D(k, ω) = 0 corresponds to the inviscid dispersion relation for gravitycapillary waves, generalized here to arbitrary depth h. This generalisation is valid in the limit of bulk-dominated dissipation for small viscosity (assumption (g)).

Ŝ(k, ω) = k N + ik • T ρ , ( 3
In the Fourier space (k x , k y , ω), the inviscid dispersion relation forms a surface with rotational invariance around the ω axis, represented by the blue green surface in figure 2. The Fourier modes (k, ω) in the immediate vicinity of the surface D = 0 correspond to propagative waves. Although the finite depth and the capillary terms are not expected to play a key role in the wrinkle generation (see section 2), we briefly recall here the main properties of the dispersion relation in the general case. Three propagation regimes may be defined: kh 1 corresponds to the shallow water regime, in which surface waves are non dispersive (finite slope at the origin in figure 2, given by the phase velocity c = ω/k = √ gh). The two other propagation regimes, lying in the deep water domain kh 1, have a phase velocity c 2 = g(1/k +k 2 c ), with c the capillary length. For k c < 1, gravity effects dominate over capillary forces, while for k c > 1 capillary forces take over, with a change of curvature in ω at k c = 1. At this inflection point the phase velocity is minimum, c min = (4gγ/ρ ) 1/4 . For the interface between air and water or the viscous aqueous solutions considered here, this minimum phase velocity is c min 22 cm/s.

In view of Eq. (3.20), a significant surface response is expected where D is small and Ŝ is large. More specifically, Eq. (3.20) predicts an excitation of the Fourier modes (k, ω) supplied by the source Ŝ that fall in the vicinity of the dispersion relation. In the simple case of pressure or stress patches rigidly travelling at a constant velocity U c in the x direction, Ŝ is non-zero for ω = k x U c , represented as a tilted plane in pink in figure 2. The intersection between this tilted plane and the surface D = 0 (shown as a black curve), which is defined provided that U c c min , is naturally where wave amplification is expected. For such idealised 'rigid' forcing, the resulting wave pattern is a collection of wakes, stationary in the frame of the source (i.e., with a convection velocity given by U c ), analogous to far-field wakes behind ships.

In the case of time-varying forcing, relevant to the problem of a turbulent wind blowing over the liquid surface, the forcing spectrum is now a continuum of Fourier modes centred around ω U c k x , in a subspace of thickness along ω given by the inverse correlation time of the fluctuations. Loosely speaking, the black line in figure 2 now has a finite thickness, allowing in principle for wave excitation even if U c < c min . We expect the wrinkles to be the unstationnary wakes generated by this turbulent forcing. For a more quantitative description, the spatio-temporal properties of the forcing Ŝ must be specified, to which section 4 is devoted.

It is worth noting that the physical picture given here for the wrinkle regime corresponds only to the linear response of the free surface, and ignores energy transfers between modes of finite amplitude and retroaction on the turbulent forcing. The regular waves, defined as the subset of Fourier modes satisfying k y = 0 (i.e., propagating in the streamwise direction), which are found experimentally at large wind velocity, probably escape from this linear description. Such regular waves first appear with a phase velocity of the order of c min , which is much smaller than the forcing velocity U c , suggesting that they originate from a retroaction of the liquid surface on the turbulent forcing. We shall return to this point in section 6.2.

Non dimensional form

We express now the surface response (3.18) in the form (2.3) inferred from the dimensional analysis. Using the boundary layer thickness δ as the characteristic length and δ/u * as the characteristic time, we introduce the dimensionless source term Ŝ † as

Ŝ † = ρ ρ a Ŝ u * δ 2 . (3.23)
The mean surface displacement ζ 2 can be related to its representation in Fourier space using Parseval's theorem,

d 2 rdt ζ 2 = 1 (2π) 3 d 2 kdω | ζ| 2 . (3.24)
Replacing ζ in Eq. (3.24) by its expression from Eq. (3.20) and Eq. (3.23) yields

ζ 2 δ 2 = ρ a ρ 2 u * 3 δ 1 (2π) 3 d 2 kdω | Ŝ † | 2 (ω 2 -ω 2 r ) 2 + ω 2 ν ω 2 , (3.25)
where ω r = g k tanh(kh) is the surface wave frequency and ω ν = 4ν k 2 is the dissipation rate. Equation (3.25) needs further specifications to be written in a dimensionless form. A choice arises for the dimensionless frequency, which can be constructed either from the characteristic time of the source, δ/u * , or from the characteristic time of the surface response, ω ν -1 . We chose the dimensionless frequencies ω = ω/ω ν , ωr = ω r /ω ν and the dimensionless wavenumber k = kδ, yielding

ζ 2 δ 2 = ρ a ρ 2 u * 3 4gν 1 (2π) 3 × d 2 kdω 1 k3 (1 + Bo -2 k2 ) tanh( kh/δ) | Ŝ † | 2 (ω -ωr ) 2 (ω/ω r + 1) 2 + (ω/ω r ) 2 .(3.26)
This equation confirms the role played by the dimensionless combination u * 3 /(gν ) introduced in section 2, and provides an analytical expression for the dimensionless function f 3 in Eq. (2.3). However, the dependency in Reynolds number is hidden in the source term Ŝ † . Hence, a quantitative description of the spectral source originating from the turbulent boundary layer in the air phase is now required.

Properties of the turbulent forcing from DNS

We compute the source term Ŝ using a set of time-resolved pressure and shear stress fields evaluated at z = 0, taken from three-dimensional DNS of a developed turbulent flow in a non deformable channel with no-slip condition at the bottom and top boundaries, and periodic boundary conditions along the streamwise and spanwise directions. In order to provide comparison with the experiments of [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF][START_REF] Paquier | Viscosity effects in wind wave generation[END_REF], which were performed in a developing boundary layer flow, we assume here that the spatio-temporal statistics of turbulence in a developing boundary layer of thickness δ at a given Re δ = δu * /ν a are equivalent to that of a channel flow of half-height H at the same value of Re τ = Hu * /ν a . Previous works have shown that this assumption is reasonably well satisfied for the flow close to the wall, z < 0.6 δ [START_REF] Jimenez | Turbulent fluctuations above the buffer layer of wall-bounded flows[END_REF][START_REF] Jimenez | Turbulent boundary layers and channels at moderate Reynolds numbers[END_REF]). In the following, we identify H to δ and we use for simplicity the same notation Re δ for the DNS and the experiments.

Boundary conditions

We first examine to what extent a canonical turbulent flow over a smooth and rigid wall with no slip boundary condition can adequately model the turbulence over a free surface. We base our analysis of the air flow on three assumptions:

• The interface is slightly deformable,

∂ t ζ| z=0 + = v z | z=0 + v x,y | z=0 + . • The interface is smooth, ζ δ ν .
• The surface drift velocity in the liquid is negligible compared to the velocity of the turbulent structures, v x | z=0 - U a . These three assumptions are corollary of the assumptions (b), (d) and (g) discussed in section 3.1. The rigid wall approximation is equivalent to the linear approximation (b). The smooth wall approximation is motivated by the typical wrinkle amplitude ζ rms , at least ten times smaller than the viscous sublayer δ ν below the wind wave threshold [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF]. It is therefore a consequence of assumption (g). Finally, the no-slip boundary condition derives from assumption (d).

Under these assumptions, we model the turbulence by a turbulent channel flow over a smooth and rigid wall with no slip boundary condition in z = 0 and z = 2H and periodic boundary conditions along x and y. The DNS configuration is sketched in figure 3, with L x , L y the streamwise and the spanwise lengths. The flow is driven by a mean streamwise pressure gradient -P a /L x . By conservation of the streamwise momentum, the mean tangential stress at the boundaries is τ a = HP a /L x , which defines the friction velocity u * = τ a /ρ a . We decompose the instantaneous pressure at the wall as the sum of a stationary pressure drop P (x) = P a (1 -x/L x ) and turbulent pressure fluctuations p(x, y, t) of zero mean. Similarly, the wall shear stress is the sum of a stationary component τ a e x and fluctuations σ(x, y, t) = σ x e x + σ y e y of zero mean. As argued before, the steady contributions are responsible for a mean flow generation in the liquid, which is not considered here. In the following, we focus on the fluctuating contributions (p, σ).

To compute the source term Ŝ = (k N + ik • T )/ρ from the (p, σ), we need to specify are the finest and coarsest spatial resolutions in the wall-normal direction in dimensionless units. ∆t + is the temporal separation between stored flow fields and Tmax is the total time simulated. Ua is the corresponding wind velocity for a boundary layer thickness δ 30 mm and a kinematic viscosity of air νa = 15 × 10 -6 m 2 s -1 (see section 2.2).

the normal and tangential stresses, where p is the air pressure and σ is the tangential wall stress in the limit of a non deformable wall. In the following, we use the conventional wall-unit notation +, {p, σ} = ρ a u * 2 {p + , σ + }, (4.6) so the dimensionless source term Ŝ † defined in Eq. (3.23) reads

N = p + 2ρ a ν a (∂ z v z ) |z=0 + (4.1) T x = σ x + ρ a ν a (∂ x v z ) |z=0 + (4.2) T y = σ y + ρ a ν a (∂ y v z ) |z=0 + . ( 4
Ŝ † (k, ω) = k p+ + i k • σ+ , (4.7)
with k = kδ.

DNS Configuration

The numerical set up configuration is sketched in figure 3, with the wall pressure and wall shear stress taken for Re δ =250. The parameters for each DNS run are summarised in table 2. The turbulent Reynolds number Re δ ranges from 100 to 550, which corresponds in the experiments of [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF] to wind speeds ranging from 1 to 5.5 m s -1 . This correspondence is obtained by equating the DNS and experimental values of Re δ , with ν a = 15 × 10 -6 m s -2 for the kinematic viscosity of air, and δ 30 mm for the local boundary-layer thickness at the x-location where measurements are carried out (see section 2.2). x as contour lines (lines are separated by increments σ + x = 1, positive for full lines and negative for dashed lines). A magnification of the snapshot by a factor 4 in each direction is also presented.

The incompressible flow is integrated in the form of evolution equations for the wallnormal vorticity and for the Laplacian of the wall-normal velocity, as in [START_REF] Kim | Turbulence statistics in fully developed channel flow at low Reynolds number[END_REF], and the spatial discretisation is desaliased Fourier series in the two wall-parallel directions and Chebychev polynomials in z. Time stepping is the third-order semi-implicit Runge-Kutta method from [START_REF] Moser | Direct numerical simulation of turbulent channel flow up to Re = 590[END_REF].

The computational box is L x ×L y ×L z = (8π, 3π, 2)δ with periodic boundary conditions along x and y directions. These spatial dimensions are larger than the minimum channel size (2π, π, 2)δ often used in numerical simulation [START_REF] Jimenez | Near wall turbulence[END_REF]. Preliminary tests showed that this large domain size is necessary to ensure the correct convergence of the Fourier integral (3.19), which is dominated by the contributions at small k.

The periodic boundary condition in time, implicitly assumed in our spectral formulation, is naturally not satisfied in the DNS data. However, spurious temporal correlations are limited by the large computational domain: the correlation time of the pressure fluctuations, of order of 20Re -1 δ δ/u * , is comfortably smaller than the transit time over the computational domain, of order of L x /U c 2δ/u * , where U c is the convection velocity (see section 4.3). The total integration time T max is chosen at least 10δ/u * (except for the largest Re δ ), to correctly resolve the lowest frequencies ω in the wave dynamics. The time step, ∆t + , is sufficiently small to resolve the fastest waves (the case Re δ = 360 has a coarser time step in order to collect statistics for a longer time period). We can therefore anticipate that the surface response will be dominated by the pressure forcing.

To compute the characteristic dimensions of the pressure and shear stress structures, we define for any field f (r, t) the spectral barycenter K,

K = K x e x + K y e y = D d 2 kdω k| f | 2 D d 2 kdω | f | 2 , ( 4
Ω = D d 2 kdω ω| f | 2 D d 2 kdω | f | 2 . (4.9)
We finally define the convection velocity U c as

U c = Ω/K x .
(4.10)

The convection velocity U c for the pressure and the shear stress, plotted in figure 5(c), lies in the range [0.5, 0.8] U a . It slightly decreases with Re δ , down to U c ≈ 0.6U a for the pressure and 0.5U a for the shear stress. These values correspond to the mean velocity at the wall-normal location z + 30 where the turbulent fluctuations are maximum [START_REF] Kim | On the structure of pressure fluctuations in simulated turbulent channel flow[END_REF]. Note however that this convection velocity is an average over Fourier components traveling at different velocities: the largest structures propagate at U c ≈ 0.8U a while the small scale structures propagate at a slightly lower value, U c ≈ 0.6U a , as observed experimentally [START_REF] Willmarth | Measurements of the fluctuating pressure at the wall beneath a thick turbulent boundary layer[END_REF][START_REF] Corcos | The structure of the turbulent pressure field in boundary-layer flows[END_REF]) and numerically [START_REF] Choi | On the space-time characteristics of wall-pressure fluctuations[END_REF].

Integrated model of the wrinkle regime

We now combine the analytical results for the surface response (section 3) with the DNS of the turbulent boundary layer (section 4) to determine the statistical properties of the wrinkles. We first compute in section 5.1 the spatio-temporal fields of synthetic wrinkles from direct integration of Eq. (3.19) using three-dimensional discrete Fourier transform, and compare them to experimental data. Although a good qualitative agreement is obtained, this direct method suffers from discretisation effects at small wavenumber. To circumvent this difficulty, we analyse the surface response in the spectral space in section 5.2, and we introduce in section 5.3 a semi-analytical method to evaluate the three-dimensional integral from its dominant contribution in the vicinity of the twodimensional resonant manifold. This method gives more insight into the physics of the wrinkles and their scaling properties.

Surface displacement computation

We first provide here a direct computation of time series of synthetic wrinkles fields from direct integration of Eq. (3.19). From the space-time Fourier transform of the wall pressure p+ (k, ω) and wall shear stress σ+ (k, ω) extracted from the DNS runs, we compute the source term Ŝ † (k, ω) from Eq. (4.7) on a discrete three-dimensional Cartesian grid (k x , k y , ω). Since the surface response occurs mainly at low wavenumber and low frequency, we perform a spectral decimation: we retain only the modes k i δ < 20 and ωδ/U a < 25.

Figure 6 shows snapshots of synthetic wrinkle fields computed using this direct method, for Re δ = 100, 250 and 550, compared to experimental measurements by [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF] for a liquid viscosity ν = 30 mm 2 s -1 and equivalent Reynolds numbers (the corresponding wind speeds are U a = 1, 2.5 and 5.5 m/s). These synthetic fields are obtained from the pressure and shear stress snapshots shown in figure 4. The wrinkles appear as disordered fluctuations, nearly isotropic at Re δ = 100, that become elongated in the streamwise direction as Re δ increases. We can note the good qualitative agreement between experimental and synthetic wrinkles for Re δ = 100 and 250 (quantitative comparisons are provided in section 5.3). By comparing figures 4 and 6, we note that the characteristic size of the wrinkles is always significantly larger than the size of the pressure and shear stress pattern from which they originate. This shift of the surface response towards larger scales is a key feature of the wrinkle regime. More precisely, the width of the wrinkles slightly decreases with Re δ , while their length remains nearly constant. The match between experimental and synthetic wrinkles is not as good at larger Reynolds number, for Re δ = 550: as the Reynolds number is increased, the DNS resolves smaller scales, while the experimental measurements reach the limit of resolution at an intermediate spatial scale, smoothing off the small scales. The experimental wrinkles at large Reynolds number appear here as a coarse-grained version of the synthetic wrinkles. Another possible source of visual discrepancy is the appearance of evanescent regular waves in the experiment, which were shown to coexist with wrinkles slightly below the wind threshold [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF]. These regular waves, which originate from an instability mechanism, cannot be captured by the present linear theory.

An interesting question is whether the wrinkles originate mostly from the pressure forcing or the shear stress forcing. The relative contribution of the two terms is illustrated in figure 7 has p + rms 4σ + rms here, see figure 5a). However, this larger amplitude is not sufficient to explain the factor 10 in amplitude between figure 7(a) and 7(b). The stronger influence of pressure also originates from the particular form of the transfer function 1/ D which tends to amplify structures of larger size. We can conclude that, although elongated in the streamwise direction, wrinkles are essentially disordered wakes generated by the nearly isotropic traveling pressure fluctuations.

It is interesting to discuss the geometry of the wrinkles in the context of the Kelvin-Mach transition observed in ship wake patterns at large Froude number [START_REF] Rabaud | Ship wakes: Kelvin or Mach angle?[END_REF][START_REF] Darmon | Kelvin wake pattern at large Froude numbers[END_REF]. We define the Froude number of a nearly isotropic pressure patch of characteristic dimension Λ x Λ y 250 δ Re -1 δ traveling at velocity U c as Fr = U c / √ gΛ. In figure 6(a,b), one has Fr 0.7, a value close to the transition Fr c 0.5 below which wakes are well described by the classical Kelvin wake pattern, of half-angle of α = sin -1 (1/3) 19.5 o . Some oblique bands can indeed be distinguished in the snapshots, reminiscent of such Kelvin wakes. In figures 6(b,c) and 6(d,e), one has Fr 2.8 and 9 respectively, for which the wakes are in the Mach-like regime, characterised by a much smaller angle α 0.2 Fr -1 ( 4.5 o and 1.5 o , respectively). The thinning of the wrinkles at increasing Re δ is therefore a signature of the decreasing angle of the wakes generated by the traveling pressure patches. We finally note that the Bond number based on the size of the pressure patches, Λ/ c 250 Bo Re -1 δ , decreases between 35 and 6 for the range of Reynolds numbers considered here. As observed in Moisy & Rabaud (2014a), the wakes are essentially in the gravity regime for these values, with weak capillary effects.

Analysis in Fourier space

An important drawback of the direct reconstruction of synthetic wrinkles from the threedimensional evaluation of Eq. (3.19) is the strong discretisation effects for small (k, ω). The most amplified Fourier modes in Eq. (3.26) occur in a narrow range of (k, ω), where the transfer function 1/ D takes large values, that is difficult to resolve numerically on a discrete grid. Although the overall shape of the wrinkle fields are robust, their amplitude shows a significant dependence on the domain size and duration of the DNS, and hence on the discretisation in (k, ω). To circumvent this discretisation issue, we perform a finer analysis of Ŝ † and ζ in Fourier space, allowing for a refined and more robust evaluation of the wrinkle properties.

Figure 8 spanwise direction, k y , and plotted in the (k x , ω) plane. By symmetry, only the two quadrants corresponding to k x > 0 are shown. The wavenumber k x is normalised by the boundary layer thickness δ, and the angular frequency is normalised by ω δ , the angular frequency of a wave of wavenumber kδ = 1. In Figure 8(a), the energy of the source Ŝ is spread over a broad band centred around the line ω = k x U c (black dashed line), where the convection velocity U c corresponds to that measured in figure 5(c). The width of the band is related to the correlation time of the turbulent fluctuations. A rigid pattern traveling at constant speed would correspond to a perfect accumulation of energy along the line ω = k x U c . The dispersion relation ω(k x ) is also plotted (red dashed line), showing that the forcing energy is mostly supplied to waves in the gravity regime (the capillary regime starts at k x δ = Bo = 14, which is outside the axis of the figure).

Representing the three-dimensional spectrum ζ(k, ω) in a two-dimensional form is delicate, because of the lack of symmetry of Eq. (3.18) in the plane (k x , k y ). We provide in figure 8(b) a two-dimensionnal representation of the spectrum, | ζ| 2 ky , using an averaging along k y as for the source Ŝ. We see that the energy of the surface response is located at smaller wavenumbers than the forcing, confirming that the wrinkles are of larger size than the pressure and shear stress structures. The resonant response is mostly contained in the ω > 0 quadrant, but it also has a significant amount of energy in the ω < 0 quadrant, indicating a small counter-propagating component. In that representation, energy accumulates around two regions: a first region surrounding the dispersion relation (red dashed line), and a second region surrounding the forcing ω = U c k x (black dashed line). This second region would suggest that a significant amount of energy is away from the resonance. However, this apparent non-resonant response is an artifact of the averaging over k y which respects the symmetry of the source but not that of the dispersion relation (see sketch of figure 2). This bias is removed in figure 9(a), showing the same spectrum now averaged in the azimuthal direction, | ζ2 | θ , as a function of k = (k 2

x + k 2 y ) 1/2 . In that representation, all the energy is now located near the dispersion relation (red dashed line). This clearly indicates that the accumulation of energy along the forcing ω = U c k x in figure 8(b) was a contribution of the Fourier components satisfying the dispersion relation with k y = 0.

The strong accumulation of energy along the dispersion relation is also present in the experiment: Figure 9(b) shows the azimutally averaged spectrum | ζ2 | θ computed from the experimental surface deformation fields for the same Reynolds number Re δ = 250. In spite of the lower spatial resolution of the experimental data, a clear accumulation of energy appears in the vicinity of the dispersion relation. We can note a slight shift of energy at frequency larger than the dispersion relation. The shift may be attributed to the surface drift current U s , which yields a Doppler-shifted dispersion relation ω g k tanh(kh) + U s k x .

The dominant response along the dispersion relation, observed both experimentally and numerically, confirms that the wrinkles are a superposition of a broad range of propagating waves. Their main specificity is their non trivial transverse structure: While regular waves correspond to k y = 0, wrinkles are characterised by wave vectors k tilted with respect to the wind direction (k y = 0), leading to elongated patterns, statistically stationary along y (because of the symmetry k y → -k y ) and propagating along x.

Integration of the resonant response

Since the main surface response occurs along the dispersion relation, we may simplify further the three-dimensional integral (3.19) by considering only the resonant response. In the limit of narrow resonance (ω ν ω r , i.e. ωr 1) and of slow varying source amplitude over the width of the resonance (

ω ν ∂ ω | Ŝ| 2 | Ŝ(k, ω r )| 2
), the integrand of Eq. (3.26) may be substituted by its second-order Taylor expansion,

| Ŝ † ( k, ω)| 2 (ω -ωr ) 2 (ω/ω r + 1) 2 + (ω/ω r ) 2 | Ŝ † ( k, ωr )| 2 1 + 2(ω -ωr ) + (4 + ω-2 r )(ω -ωr ) 2
(5.1)

+ O((ω -ωr ) -2 ).
Note that the odd term (ωωr ) of the expansion does not contribute to the wrinkle amplitude as it cancels out upon integration along ω.

The validity of this approximation is illustrated in the inset of figure 9(a): the compen-

sated spectrum 2πk | ζ(k, ω)| 2
θ , plotted as a function of the angular frequency ω for a fixed wavenumber, kδ = 2.5 (white curve), is indeed very close to its second order Taylor expansion (red dashed line). The limit of narrow resonance breaks down for ω ν ≈ ω r , which would occur around kδ = 10 for viscosity ν = 30ν water . In practice, most of the energy is contained in the region kδ < 10, so that the full integral of Eq. (3.26) can be safely approximated by the limit of narrow resonance.

In that limit, the integral along ω in Eq. (3.26) can be performed analytically, and we obtain

| ζ2 | ω (k) = dω| ζ(k, ω)| 2 = ρ a ρ 2 u * 3 16gν 1 (2π) 2 W ( k) | Ŝ † ( k, ωr )| 2 , (5.2)
where we introduce the weighting factor

W ( k) = 1 k3 (1 + Bo -2 k2 ) tanh( kh/δ) . (5.3)
This weighting factor originates from the expression of the dispersion relation, and is responsible for the shift towards larger scales of the surface response. This wavelength shift depends on the values of the Bond number Bo and the dimensionless depth h/δ. For k 1/h, the waves are in a shallow water regime, and the weighting factor is W ( k) δ/(h k4 ). For 1/h k 1/ c , the waves are in the gravity regime and W ( k) 1/ k3 . Finally, for k 1/ c , the waves are in the capillary regime and W ( k) Bo 2 / k5 . In all cases, this weighting factor tends to enhance the low-k content of the forcing. In practice, most of the energy is contained in the deep-water gravity regime, so that W ( k) 1/ k3 is the most relevant weighting factor in the wrinkle problem.

From Eq. (5.2), the surface response can be computed using a linear interpolation of Ŝ † (k, ω) on the dispersion relation manifold D(k, ω) = 0. The interpolation is performed along the direction ω using 20 mesh points in the range [0.8ω r , 1.2ω r ] and the same mesh size in k than the original mesh size, and the value at the resonance is taken as the average over the 20 computed values of ω for each k. This procedure reduces the computational cost by a factor up to 100 for the highest Reynolds number. Its main advantage is to remove the evaluation of the full three-dimensional integral which presents a highly peaked resonance near ω = ω r . As a consequence, the properties of the wrinkles are found to be less sensitive to the spatial and temporal discretisation of the DNS data.

Figure 10 shows the Fourier spectrum of the surface deformation (5.2) for increasing values of the Reynolds number Re δ . The energy is essentially contained in a zone located at relatively small k, with a symmetric distribution along k y = 0. It essentially follows the intersection of the dispersion relation manifold and the forcing plane ω = k x U c , given by

gk tanh(kh)(1 + (k c ) 2 ) = (k x U c ) 2 ,
(5.4) plotted in red line using for each Re δ the convection velocity measured in figure 5(c). While the energy remains mostly located in the range k x δ ∈ [0, 2], it significantly broadens in the k y direction as Re δ increases, i.e. the angle of the phase velocity gradually departs from the direction of the wind. This increase of k y /k x is the spectral signature of the thinning of the wrinkles in the spanwise direction. More specifically, we can compute the spectral barycenter K of | ζ| 2 from Eq. (4.8), shown as circles in figure 10. The corresponding streamwise and spanwise sizes, defined as Λ x = 2π/K x and Λ y = 2π/K y , are plotted in figure 11(a,b), and compared to the experimental data of [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF].

The almost constant streamwise size Λ x /δ = 6.7 ± 0.7 and the decreasing spanwise size Λ y /δ are qualitatively recovered. It is worth noting that, in spite of their large streamwise extent compared to the liquid depth, wrinkles are essentially deep-water waves: their wavenumber k is dominated by the spanwise component k y , for which we have tanh(kh) > 0.97 in this range of Re δ . Last but not least, the computed convection velocity, shown in figure 11(c), closely follows the experimental data: at small Re δ , the normalized convection velocity U c /U a of the wrinkles is close to that of the pressure forcing (U c /U a 0.7, blue squares), but it decreases significantly as Re δ increases. This decreasing convection velocity for the wrinkle is a consequence of the propagation angle of the Fourier modes composing the wrinkles, which gradually departs from the direction of the wind as Re δ increases.

The correct agreement between synthetic and experimental wrinkles can be extended to the range of viscosity ν water < ν < 10 3 ν water , since both the experimental measurement and the theory do not exhibit variation in viscosity on sizes and convection speed. At large Reynolds number, Re δ > 400, the sharp decrease in Λ x and increase in Λ y found experimentally is not reproduced by the computation. This sharp evolution is associated to the instability that gives rise to the regular waves, which cannot be captured by the present linear model.

We finally turn to the scaling of the wrinkle amplitude as a function of the friction velocity u * and liquid viscosity ν . By considering only the resonant contribution, the wrinkle amplitude ζ rms is obtained from Eq. (5.2) by summing | ζ2 | ω over k,

ζ 2 rms = ρ a ρ 2 u * 3 16gν 1 (2π) 2 d 2 k W ( k) | Ŝ † ( k, ωr )| 2 .
(5.5)

Figure 12(a) shows the wrinkle amplitude ζ rms /δ as a function of Re δ . The experimental scaling in ζ rms ∝ u * 3/2 is well reproduced by the the synthetic wrinkles, but with an amplitude twice larger. This discrepancy probably originates from the high sensitivity to the low-wavenumber content of the forcing. The contribution of the largest scales to the wrinkle amplitude ζ rms can be different experimentally and numerically for two reasons. First, the measurements were carried out on a window smaller than the channel width, so that the smallest wavenumbers may be poorly estimated. Second, the numerical simulations are performed in a box with periodic boundary conditions, so the largest scales could be different from that of a true developing turbulent boundary layer. Finally, we show in figure 12(b) the experimental wrinkle amplitude ζ rms as a function of the dimensionless liquid viscosity ν g/u * 3 for a fixed value of u * corresponding to Re δ = 180. The data are in good agreement with the analytical prediction ζ rms /δ ∝ (ν g/u * 3 ) -1/2 . We can conclude that the dimensionless function f 4 introduced in Sec. 2.1 is essentially independent of the Reynolds number in the range Re δ ∈ [100, 550]. The additional dependencies in Bo and h/δ, explicitly considered in the derivation, cannot be tested against experiments, which were performed for fixed surface tension and liquid depth (Bo 14 and δ/h 1.2, see section 2.2). Ignoring these dependencies in Bo and h/δ, the function f 4 reduces to a constant, and Eq. (2.4) simply writes:

ζ rms δ C ρ a ρ u * 3/2 (gν ) 1/2 .
(5.6)

The numerical factor, fitted in figure 12(b), is C 0.022 ± 0.005.

6. Connection with the inviscid resonant theory of [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF] The present theory focuses on the statistically steady wrinkle regime, of amplitude governed by the liquid viscosity. This regime is the asymptotic state of the surface deformation, reached when the energy input by the turbulence forcing is balanced by the viscous dissipation. Before this energy balance is reached, a transient growth regime must take place, where viscous dissipation can be neglected. We show here how this inviscid growth regime, previously investigated by [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF], naturally asymptotes towards the viscous saturated wrinkle regime described in the present theory, provided that the wrinkles remain of small amplitude.

Temporal growth of wrinkles

In the previous sections, all fields were assumed statistically homogeneous and stationary, allowing for a space-time Fourier description. To describe the transient growth that precedes this steady regime, we still assume here homogeneity but we relax the stationary assumption. Only a spatial Fourier transform of the dynamical equation is then performed: (6.1b) noting ζ the spatial Fourier transform for ζ and similarly for the pressure and stress fields. The dynamics is now governed by a Langevin equation [START_REF] Langevin | Sur la théorie du mouvement brownien[END_REF][START_REF] Pottier | Non Equilibrium Statistical Physics[END_REF] for the stochastic wave amplitude ζ, which can be derived following an approach similar to the derivation of section 3,

ζ(k, t) = d 2 rζ(r, t)e -ik•r (6.1a) ζ(r, t) = (2π) -2 d 2 k ζ(k, t)e ik•r ,
∂ tt ζ(k, t) + 4ν k 2 ∂ t ζ(k, t) + g k ζ(k, t) = -k p0 (k, t) -ik • σ0 (k, t). ( 6
c (k) ∼ (k • U k ) -1
, characterising the temporal correlation of the turbulent structures, and a slow dissipation time τ ν (k) ∼ (ν k 2 ) -1 , associated to viscous dissipation in the liquid. We assume here for simplicity that the convection velocity U k is the same for all k, and given by the global convection velocity U c e x . In the intermediate-time regime τ c t τ ν , each mode k corresponds to an essentially undamped oscillator forced by an uncorrelated noise, resulting in a linear growth of | ζ(k, t)| 2 , with a k-dependent growth rate governed by the corresponding Fourier component of the pressure forcing (the shear stress forcing may be ignored during this quasi-inviscid growth). The mean square wave amplitude, integrated over all modes, similarly grows linearly in time, resulting in classical result (1.2) of [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF]. For time larger than the slowest (largest-scale) growing mode τ ν δ 2 /ν , all Fourier components are saturated, and the asymptotic mean square amplitude can be simply estimated by setting t τ ν in Eq.(1.2): with p 2 ∝ ρ 2 a u * 4 and U c ∝ u * , we recover the result of Eq. (5.6). The wrinkle regime described in this paper therefore naturally arises as the viscous-saturated asymptotics of the inviscid growth theory of [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF].

Onset of regular waves

A key assumption in our theory, as well as in the inviscid resonant theory of [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF], is the absence of feedback of the wrinkles deformations on the turbulence in the air. In other words, the growth and saturated regimes sketched in figure 13 hold only provided that the wrinkle amplitude remains small compared to the thickness of the viscous sublayer δ ν in the turbulent air flow: this is assumption (g), used to derive the linearised surface response. The questions that naturally arise now are what is the maximum wrinkle amplitude before the breakdown of this assumption, and whether this breakdown could be related to the onset of regular waves.

The breakdown of the decoupled dynamics hypothesis (assumption (g)) can be expected when the amplitude of the wrinkles reaches a fraction of the viscous sublayer thickness, ζ Aδ ν , with δ ν = ν a /u * and A a numerical factor. Using Eq. (5.6), this criterion is satisfied for a friction velocity u * beyond a critical value, law (6.4) yields A 0.11 ± 0.02. This good match with the scaling ν 1/5 suggests that regular waves could be triggered by an instability originating from the feedback of the wrinkles on the air turbulence: once the wrinkle amplitude reaches Aδ ν , the pressure and shear stress fluctuations in the boundary layer are no longer that of a no-slip flat surface, but acquire a spatio-temporal structure reflecting the shape of the surface. In turn, this spatio-temporal phase coherence between the wave field and the forcing could enhance the energy transfer, leading to the exponentially growing waves found in experiments. In this scenario, wrinkles appear as the natural base state from which regular waves grow. This scenario cannot be tested by the present theory, which ignores such coupling between the liquid and the air phases.

u * c = A C 2/5 ρ ρ a 2/5 gν ν 2 a δ 2 1/5 . ( 6
It may be noted that this tentative criterion for wave onset suggests that the turbulent boundary layer becomes sensitive to the surface roughness for rms amplitude of order of 0.1δ ν . Such roughness is surprisingly small: the peak of turbulent kinetic energy in a boundary layer is at 15δ ν , and the boundary layer is essentially a laminar shear flow up to 4δ ν . Boundary layer turbulence over a wavy no-slip wall is indeed essentially unaffected by rigid wall roughness up to 4δ ν [START_REF] Schlichting | Boundary Layer Theory[END_REF][START_REF] Jimenez | The large-scale dynamics of near-wall turbulence[END_REF]). The relatively small wrinkle amplitude found here for the growth of regular waves probably originates from the specific phase coherence of the surface waves and the pressure perturbations they induce: this phase coherence possibly enables an optimal energy transfer, and hence an exponential growth of regular waves even from very fine seeding wrinkles.

Conclusion

In this paper a spectral theory is derived to describe the surface deformations of small amplitude under arbitrary normal and tangential stresses applied at the air-liquid interface (wrinkle regime), assuming no feedback of such deformations on the air flow.

The key result of the paper is the demonstration of the scaling for the wrinkle amplitude, ζ/δ (ρ a /ρ )u * 3/2 /(gν ) 1/2 , in good agreement with the experimental findings of [START_REF] Paquier | Surface deformations and wave generation by wind blowing over a viscous liquid[END_REF][START_REF] Paquier | Viscosity effects in wind wave generation[END_REF]. This theory corresponds to the viscous-limited asymptotic steady state of the inviscid resonant mechanism proposed by [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF], and provides an appropriate description of the surface deformations for wind velocity below the onset of regular waves.

A significant improvement of the present theory is the quantitative description of the fraction of energy supplied by the pressure fluctuations that is located near the resonance. As already pointed out by [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF], only the pressure fluctuations of space-time correlations matching the dispersion relation contribute to the surface deformations. Detailed knowledge of the space-time Fourier spectrum of the pressure and shear stress fluctuations in a turbulent boundary layer, which was not available at the time of [START_REF] Phillips | On the generation of waves by turbulent wind[END_REF], was used here to close the problem by determining numerically the dependence of the wrinkle amplitude with the governing parameters. The wrinkle regime therefore provides an interesting configuration where the effect of a turbulent forcing on a dispersive wave system can be exactly computed. A similar approach was recently proposed for waves generated on a viscoelastic compliant coating [START_REF] Benschop | Deformation of a linear viscoelastic compliant coating in a turbulent flow[END_REF].

We have shown that the wrinkles below the wave onset correspond to a superposition of uncoherent wakes mostly originating from the pressure fluctuations traveling in the turbulent boundary layer (the shear stress fluctuations are found to provide a negligible contribution to the wrinkles). The thinning of the wrinkles in the spanwise direction as the wind velocity increases is reminiscent of the decrease of the wake angle for a finitesize moving disturbance found in the classical ship wake problem [START_REF] Rabaud | Ship wakes: Kelvin or Mach angle?[END_REF][START_REF] Darmon | Kelvin wake pattern at large Froude numbers[END_REF]Moisy & Rabaud 2014b). This mechanism could be related to the surprisingly large cross-wind wave slopes found in ocean observations [START_REF] Munk | An inconvenient sea truth: Spread, steepness, and skewness of surface slopes[END_REF].

A remarkable property of the wrinkle is that their characteristic size are governed by the (outer) boundary layer thickness δ, although they originate from pressure patches of characteristic size governed by the (inner) viscous sublayer thickness, δ ν δRe -1 δ . This is because the liquid surface response integrates the pressure forcing, resulting in a systematic shift towards the upper bound of the energy-containing range [δ ν , δ] of the forcing. The wrinkles are therefore essentially governed by the largest scales of the turbulent flow. As a consequence, the detailed statistics of the wrinkles is expected to depend on the geometry of the forcing, making a fine comparison between simulations, laboratory and outdoors experiments difficult.

The present theory neglects the influence of surface drift current. Although significant effects of the liquid current are certainly present for the onset and amplification of regular waves, we expect weak influence of drift current on the wrinkle regime. Indeed, the dominant Fourier modes are oriented along the wind direction for regular waves, while they are nearly normal to the wind in the wrinkle regime. As a consequence, only a small Doppler shift will arise in the wrinkle regime, even in the presence of a strong drift current along the wind direction. We therefore expect robust properties of the wrinkles, that may be extended to liquids of small viscosity and large depth, provided that the flow in the liquid remains laminar. As the laminar condition is not satisfied for the air-sea interface, the influence of water current on the wrinkles in oceanographic conditions still deserves further analysis In spite of this limitation, the implications of the present work for physical oceanography are important. In particular, the transition towards regular (quasi-monochromatic) waves as the wind velocity is increased raises the question of the role of the wrinkles as a base state for wave amplification. The experiments of [START_REF] Paquier | Viscosity effects in wind wave generation[END_REF] suggest that the regular waves are triggered when the wrinkle amplitude reaches a fraction of the viscous sublayer thickness. Beyond that amplitude, the feedback of the surface roughness on the turbulent boundary layer can no longer be neglected. This provides a criterion for the wave onset, u * c ∝ ν 1/5 , which is consistent with experiments performed in viscous liquids, and for which no explanation has been proposed so far. Describing the surface deformations for a wind velocity above this threshold is beyond the scope of the present linear theory, in which such coupling between the liquid and the air phases is ignored. This deserves further investigation, as it could renew our understanding on the onset of wave generation.

Using the relation m 2 = k 2 -iω/ν , we obtain

-g k ζ + m 2 + k 2 m 2 -k 2 k ρ p0 + 2mk m 2 -k 2 ν Bz = k N ρ (A 6) 2k 2 m 2 -k 2 k ρ p0 + m 2 + k 2 m 2 -k 2 ν Bz = ik • T ρ . ( A 7)
This expression can be simplified by replacing p0 by its expression from Eq. (A 4)

(ω 2 -g k) ζ - 2m m + k ν Bz = k N ρ - ik • T ρ (A 8) ω 2 ζ + m 2 -k 2 2k 2 ν Bz = m 2 -k 2 2k 2 ik • T ρ . ( A
(m + k)(m 2 -k 2 ) ω 2 ζ + (ω 2 -g k) ζ = k N ρ + m -k m + k ik • T ρ . ( A 11)
Factorizing by ζ, we obtain the final expression

ω 2 -g k + 4iν ωk 2 + 4ν 2 k 3 (m -k) ζ = k N ρ + m -k m + k ik • T ρ , ( A 12)
which is Eq. (3.17).

  and using Eqs. (3.10)-(3.11), we obtain the expression of the velocity

  .15) where g = g + γk 2 /ρ is the modified gravity and Bz = (κ × Ω0 (k, ω)) • e z satisfies ∆v z = -B z . The component B z is associated to the viscous dissipation of the vertical component v z of the velocity field. Using the kinematic condition ∂ t ζ = (v z ) z=ζ in the small perturbation limit, we obtain the relation between ζ, p0 and Bz , Bz can be expressed using Eqs (3.14) and (3.15). From Eqs. (3.14), (3.15) and (3.16) we obtain by algebric manipulations (see details in Appendix A) the following expression

  .18) Equation (3.18) is the corner stone of this paper: it relates the Fourier component of the displacement field ζ(k, ω) to the Fourier components of the applied normal and tangential stresses N and T . The surface response in the physical space is finally obtained by applying the inverse Fourier transform (3.5b) to Eq. (3.18),

  Figure 2. Representation in the Fourier space of the dispersion relation ω = g k tanh(kh) of gravity-capillary waves (blue-green surface) and the forcing (pink plane). The forcing here corresponds to a source traveling at constant velocity in the x direction, ω = Uckx. The intersection between the two surfaces (black line) is where the maximum wave amplitude is expected.

  .21) and the (inverse) spectral convolution kernel D as D(k, ω) = ω 2 -g k tanh(kh) + 4iν ωk 2 .(3.22)

  .3)The relative importance of the viscous contributions depends on the surface deformation and the magnitude of the turbulent air flow. For a surface displacement ζ of characteristic wavenumber k and convection speed U c , spatial gradients of v z scale as ∂ x,y,z v z ≈ ζU c k 2 . For a turbulent boundary layer of friction velocity u * , the vertical gradient of horizontal velocity scales as ∂ z v x,y ≈ u * /δ ν . The ratio ∂ x,y v z /∂ z v x,y is thus given by ζδ ν u * k 2 /U c , where U c ≈ 12u * . If we consider a wave of wavelength Λ, using assumption (b), ζ Λ, and assumption (g), δ ν Λ, we have in practice ∂ x,y,z v z ∂ z v x,y in the wrinkle regime. The expressions of N and T thus reduce to

Figure 3 .

 3 Figure 3. Sketch of the numerical setup. The pressure and the shear stress on the wall plane z = 0 taken from Direct Numerical Simulation of a turbulent channel flow are applied to the surface of a viscous liquid. The DNS is performed on a domain (Lx, Ly, Lz) = (8π, 3π, 2)δ. The snapshot illustrates the turbulent field in z = 0 for the case Re δ = 250. Pressure p + = p/(ρau * 2 ) is shown in color, and longitudinal shear stress σ +x as contour lines (lines are separated by increments σ + x = 1, positive for full lines and negative for dashed lines). A magnification of the snapshot by a factor 4 in each direction is also presented.

Figure 4 .

 4 Figure 4. Snapshots of the pressure and shear stress fields at the surface, at three values of Re δ . Only a subdomain [0, 12δ] × [0, 8δ] is shown. The pressure p + = p/(ρau * 2 ) is shown in colour, and the longitudinal shear stress σ +x as contour lines, such that |σ + x | = 0.5i with integer i (positive for full lines and negative for dashed lines).

4. 3 .

 3 Figure4shows snapshots of the instantaneous wall pressure p + and streamwise wall shear stress σ +x for Reynolds number Re δ = 100, 250 and 550. Increasing the Reynolds number naturally decreases the size of the structures. These snapshots show that the shear stress patterns tend to be elongated in the streamwise direction, whereas the pressure patterns are nearly isotropic in the (x, y) plane. These elongated shear stress patterns are a classic signature of the streamwise streaks in the near-wall region of the boundary layer, whereas the nearly isotropic pressure patterns are related to the imprint created at the wall by the cores of the vortices[START_REF] Jimenez | Near wall turbulence[END_REF].The intensities and sizes of the pressure and shear stress fluctuations are quantified in figure5as a function of the Reynolds number. The pressure r.m.s. (figure5(a)) is typically 2 to 5 times larger than the shear stress r.m.s. and increases weakly with Re δ , while the shear stress r.m.s. remains nearly constant over the range of Re δ considered.

Figure 5 .

 5 Figure 5. Statistics of pressure and shear stress fluctuations as a function of Re δ . (a) Root-mean-squared pressure p + rms (•) and stress σ + rms = (σ + x rms 2 + σ + y rms 2 ) 1/2 ( ). (b) Characteristic streamwise and spanwise lengths Λx/δ ( , ) and Λy/δ ( , ) of pressure and shear stress, computed from the spectral barycenters (4.8); for σ, the values of Λ are averaged over the two components σx and σy. Structure sizes decrease as Λ ∝ Re -1 δ (black solid line). (c) Mean convection velocity Uc/Ua for pressure (•) and shear stress ( ), computed using Eq. (4.10).

  .8) where f (k, ω) is the Fourier transform of f , and D = {(k x , k y , ω)|k x > 0, k y > 0} is the domain of integration. The mean structure size in the streamwise and spanwise directions, defined as Λ x = 2π/K x and Λ y = 2π/K y , are plotted in figure5(b) for the pressure and the shear stress. Λ x and Λ y both decrease as Re δ -1 , indicating that they scale as the (inner) viscous sublayer thickness δ ν . The sizes normalized by δ ν are Λ +x shear stress patches. Similarly, we can define the frequency barycenter Ω of a field f as

  , obtained at the intermediate Reynolds number Re δ = 250. Figure 7(a) and (b) show snapshots of surface deformation computed using the pressure contribution only ( Ŝ † = k p+ ) and shear stress contribution only ( Ŝ † = i k • σ+ ), while figure 7(c) combines the two contributions. Pressure clearly dominates the wrinkle generation. This important result can be primarily attributed to the larger rms amplitude of pressure (one Re δ = 250, exp. (d) Re δ = 250, num. (a) Re δ = 100, exp. (b) Re δ = 100, num.

Figure 6 .

 6 Figure 6. Comparison between experimental surface displacement field ζ measured by Paquier et al. (2015) (left) and synthetic wrinkle fields computed from DNS data with Eq. (3.19) (right). Only a subdomain of the synthetic wrinkle field is shown, to match the size of the experimental domain. Experimental data (left) for Ua = 1 m/s (a), 2.5 m/s (c), and 5.5 m/s (e), which corresponds to Re δ = 100 (b), Re δ = 250 (d), Re δ =550 (f).

Figure 7 .

 7 Figure 7. Analysis of the respective role of pressure and shear stress fluctuations on the wrinkle generation. Synthetic wrinkle field snapshot obtained for Re δ = 250 from (a) the pressure contribution p only; (b) the shear stress contribution σ only; (c) both pressure and shear stress. The full field is almost indistinguishable from the field (a), showing that the main contribution originates from pressure fluctuations.

Figure 8 .

 8 Figure 8. (a) Space-time spectrum | Ŝ| 2 ky of the source term for Re δ = 250. (b) Space-time spectrum of surface displacement | ζ(ω, kx)| 2 ky , computed from Eq. (3.18) for a liquid viscosity corresponding to ν = 30 mm 2 /s. Black dashed line: ω = kxUc, where Uc is the convection velocity. Red dashed line: dispersion relation. kx is normalised by the boundary layer thickness δ, and ω is normalised by the frequency ω δ corresponding to a wave of wavenumber kδ = 1.

Figure 9 .

 9 Figure 9. (a) Space-time spectrum of surface displacement | ζ| 2 θ averaged over the azimuthal direction for ν = 30 mm 2 /s, computed from Eq. (3.18) with the DNS data at Re δ = 250. The three dimensional spectrum is averaged along the direction θ to respect the symmetry of the dispersion relation. The black line shows the forcing ω = Uck, and the red dashed line is the dispersion relation. Inset: Spectrum for a given wavenumber kδ = 2.5 computed from DNS (solid white line), and its Taylor expansion (5.2) around the corresponding resonant frequency ωr (red dashed line). The peak half-width is given by the viscous frequency, ων = 4ν k 2 0.03ω δ . (b) Space-time spectrum | ζ| 2 θ computed from the experimental data from Paquier et al. (2015), for Ua = 2.5 m/s and ν = 30 mm 2 /s.

Figure 10 .

 10 Figure 10. Spectrum of the surface displacement | ζ| 2 integrated along the dispersion relation for Re δ = 100, 180, 250, 360 and 550. The resonant curve, defined as the intersection between the forcing plane ω = kxUc and the surface ωr = g k tanh(kh), is shown in red line. The spectral barycenter (Kx, Ky) is shown (•) for each quadrant kx > 0.

Figure 11 .

 11 Figure 11. Sizes and convection velocity of the wrinkles as a function of Re δ . Comparison between the synthetic wrinkles (red filled symbols) and the experimental wrinkles of Paquier et al. (2015) (black open symbols). (a) Streamwise size Λx/δ, (b) spanwise size Λy/δ, and (c) convection velocity Uc/Ua computed from Eqs. (4.8), (4.9) and (4.10). The blue squares show the convection velocity of the pressure forcing (see figure 5(c)).

Figure 12 .

 12 Figure 12. Scaling of the wrinkle amplitude ζrms/δ with respect the Reynolds number Re δ and the liquid viscosity ν . (•) experimental data; (•) synthetic wrinkles. (a) ζrms/δ as a function of Re δ for ν = 30 mm 2 /s. (b) ζrms/δ as a function of the dimensionless number ν g/u * 3 for Re δ = 180 and liquid viscosity in the range 1 -600 mm 2 s -1 . A best fit with the analytical prediction (5.6) is shown in black solid line, yielding C 0.022.

Figure 13 .

 13 Figure 13. Sketch of the surface deformation amplitude ζ/δ as a function of time for a given Re δ and two liquid viscosities ν . Starting from a purely flat interface, three regimes follow one another: linear response at short time, for t τc, where τc δv/u * is the correlation time of the turbulent structures; Phillips (1957) regime of quasi-inviscid resonant growth at intermediate time, for τc t τν , where τν ν /δ 2 is the viscous timescale; asymptotic viscous-saturated wrinkle regime at long time for t τν . If the wrinkle amplitude reaches a given fraction of the viscous sublayer thickness (star symbol, for ζ/δ Re -1 δ ) before the viscous saturation, the linear assumption of the present theory breaks down, resulting in the possible triggering of a wave instability.

  9) Multiplying Eq. (A 9) by (mk)m/k gives 4mk 2 (m + k)(m 2 -k 2 )

Table 2 .

 2 Re δ ∆x + ∆y + ∆z + min ∆z + max ∆t + Tmaxu * /δ Ua (m s Details of the DNS turbulent channel flow. ∆x + and ∆y + are the spatial resolutions in terms of Fourier modes before dealiasing. ∆z + min and ∆z + max

	-1 )

  .2) Each Fourier component k describes a linear damped oscillator forced by a stochastic noise given by the corresponding Fourier component of the applied pressure and shear stress fields. Such Langevin equation with short-time temporal correlations in the noise term exhibits three regimes, sketched in figure13: ballistic motion at short time ( ζ ∝ t), diffusive process at intermediate time ( ζ ∝ t 1/2 ), and asymptotic regime governed by cumulative effect of viscosity at large time ( ζ ∝ t 0 ). The intermediate-time regime, defined only for liquids of small viscosity, corresponds to the inviscid resonant theory of[START_REF] Phillips | On the generation of waves by turbulent wind[END_REF], whereas the large-time saturated regime corresponds to the wrinkles.More specifically, we can introduce for each Fourier component k a fast correlation time τ

Acknowledgements

The authors thank L. Deike, C. Garrett, J. Jiménez, W. Munk, C. Nové-Josserand, A. Paquier and E. Raphaël for fruitful discussions. This work was supported by the project "ViscousWindWaves" (ANR-18-CE30-0003) of the French National Research Agency, and by the LabeX LaSIPS (ANR-10-LABX-0040-LaSIPS) managed by the French National Research Agency under the "Investissements d'avenir" program (ANR-11-IDEX-0003-02). A.L.D. acknowledges the support from the Office of Naval Research under Grant #N00014-16-S-BA10.

Appendix A. Detailed calculation of Eq. (3.17)

We detail here the calculation steps of section 3.3 which establishes the final expression of Eq. (3.17). After introducing the Fourier transforms of p , Ω and v and using the boundary conditions in z = 0 we obtain the system of equation

where g = g + γk 2 /ρ is the modified gravity and Bz = (κ × Ω0 (k, ω)) • e z is the non potential flow part of vz , which satisfies ∆v z = -B z . In Fourier space, we can evaluate F{(∂ zz v z ) z=0 } by derivating twice the expression of vz with respect to z,

Using the kinematic condition ∂ t ζ = (v z ) z=ζ in the small perturbation limit, we obtain the relation between ζ, p0 and Bz ,

Replacing F{(∂ zz v z ) z=0 } by its expression in Eq. (3.15) yields 5)