Table Of Content

Freely decaying weak turbulence for sea surface gravity waves M. Onorato, A. R. Osborne, M. Serio Dip. di Fisica Generale, Università di Torino, Via Pietro Giuria 1, 10125 Torino, Italy D. Resio Coastal and Hydraulics Laboratory, U.S. Army Engineer Research and Development Center, Halls Ferry Road, Vicksburg, MS 39180, USA A. Pushkarev Waves and Solitons LLC, 918 W. Windsong Dr., Phoenix, AZ 85045, USA 2 V. Zakharov 0 0 Landau Institute for Theoretical Physics, Moscow, 117334, Russia and Department of Mathematics, University of Arizona, 2 Tucson, AZ 85721, USA n C. Brandini a J La.M.M.A., Regione Toscana, Via A. Einstein 35/b 50013 Campi Bisenzio - Firenze, Italy 1 Westudy numerically thegeneration of power laws in theframework of weak turbulencetheory 1 forsurfacegravitywavesindeepwater. Startingfromarandomwavefield,weletthesystemevolve numericallyaccordingtothenonlinearEulerequationsforgravitywavesininfinitelydeepwater. In ] D agreement with the theory of Zakharov and Filonenko, we find the formation of a power spectrum characterized by a power law of theform of |k|−2.5. C . n After the pioneering work by Kolmogorov [1] on the small-amplitude waves. Thus, contrary to fully devel- li equilibrium range in the spectrum of an homogeneous oped turbulence, it leads to explicit analytical solutions n and isotropic turbulent flow, there have been a number provided some assumptions are made. The first exper- [ of studies on cascade processes in many other fields of imental support of the theory for surface gravity waves 1 classical physics such as plasma physics, magnetohydro- was made by Toba [6] who was completely unaware of v dynamicsandoceanwaves. Forsurfacegravitywavesthe the paper by Zakharov and Filonenko. He reformulated 7 first seminal theoretical work was done by O.M. Phillips the Phillips’ equilibrium range law in the following way: 1 in 1958, [2]. Using dimensional arguments, he argued F(ω)=βgu∗ω−4, where β should now be a universaldi- 0 that the frequency spectrum in the inertial range was of mensionlessconstant. AftertheworkbyToba,successive 1 the form F(ω) = αg2ω−5, where α was supposed to be experimentalobservationoftheω−4 lawhavebeenmade 0 2 an absolute constant and g is gravity. Even though in by a number of authors, see for example [7–10]. 0 the introduction of Phillips’ paper it was stated that “a Eventhoughthereisaconsensusonthisresult,itmust / necessary condition for the equilibrium range over a cer- n be stressed that so far the verification of the theory has tain part of the spectrum is the appreciable non-linear i never been established from first principles and more- nl interactions among these wave-numbers” (from [2]), his over the mechanisms that lead to the power law ω−4 : arguments were based on the geometricalfeatures of the are not universally recognized: geometrical aspects re- v free surface elevation. One of his basic assumptions was i lated to wave breaking, without invoking the nonlinear X that the only variable of interest was gravity, while the wave-wave interaction mechanism, are still retained by r friction velocity, u∗, was not supposed to be involved in many oceanographers as fundamental for generating an a thespectralrelation,limitingthepossibilityforacorrect ω−4 powerlaw. ConfirmationoftheZakharov-Filonenko dimensional analysis. solution to the kinetic equation has being given through numerical simulations of the kinetic wave equation itself Some years later Zakharov and Filonenko [3] estab- [11] [12], solving exactly the so called S term. Never- lished that in infinite water the direct cascade should nl theless,itmustbeunderlinedthatthekineticequationis produce apowerspectrumofthe surfaceelevationofthe form P(|k|) ∼ |k|−2.5 that corresponds, using the linear derived from the primitive equations of motion under a dispersionrelationininfinite depth, toanω−4 frequency number ofhypotheses (see for example [13]), therefore it cannotbeconcludedapriorithatpowerlawsolutionsof power spectrum: the result was found as an exact solu- thekineticequationarealsosharedbythefullynonlinear tionofthekineticwaveequation(see[4]). Thetheoryde- wave equations. velopedisknownas”weak”or”waveturbulence”andhas manyimportantapplicationsindifferentfieldsofphysics One way to verify weak turbulence theory is to per- suchashydrodynamics,plasmaphysics,nonlinearoptics, form direct numerical simulations of the primitive equa- solid state physics, etc. see [5]. It is called weak turbu- tions of motion. The numerical confirmation of the the- lence because it deals with resonant interactions among ory for gravity waves propagating on a surface has not 1 been an easy task (for capillary waves see [14], for one is of infinite depth): dimensional wave turbulence see [15,16]), basically because of the intrinsic difficulties of the computation of 1 2 2 2 2 2 ψ +gη+ ψ +ψ −(φ | ) (1+η +η ) =0 (1) theboundaryconditions. Differentnumericalapproaches t 2(cid:20) x y z η x y (cid:21) have been used for integrating the fully nonlinear surface gravity waves equations (see [17] for a review). The 2 2 numerical methods based on volume formulations show η +ψ η +ψ η −φ | (1+η +η )=0, (2) t x x y y z η x y very interesting results, in particular they are capable of modeling in a quite appropriate way wave breaking. Themajordifficultyfornumericalsimulationsofthesys- Unfortunately they have the disadvantage that they re- tem (1)-(2) consists in that we have to compute the quire large computational resources, and therefore are derivatives of φ with respect to z on the surface η. This not suitable for long time numerical simulations. For ir- problem can be overcome if we express the velocity po- rotational and inviscid flows boundary formulations are tential ψ(x,y,t) as a Taylor expansion around z = 0. usually preferred: only the surface is discretized reduc- Inverting asymptotically the expansion one can express ing the dimension of computation (from three to two). φz|η as an expansion of derivatives of ψ(x,y,t) that can The Higher-Order Spectral Methods (HOS), indeed the thenbecomputedusingtheFastFourierTransform,sim- methodusedinournumericalsimulations,introducedin- plifying notably the computation. This is nothing other dependently by West et al. [18] and by Dommermuth et than a different way for formulating the HOS method. al. [19], belongs to this second approach(see also the re- We underline that this is the same approach that has centworkbyTanaka[20]). Veryrecentlythreenewmeth- originally been used in [4] for deriving analytically the ods havebeenproposedas verypromisingfor simulating equation that is usually known as the “Zakharov equa- waterwaves[21–23]. Resultsusingthesenewapproaches tion”. The order of the simulation can be decided a pri- on turbulent cascades are still to be completed. ori and depends on how many terms are retained in the In this Letter we establish numerically, using a HOS Taylor expansions; in our numerical simulations we con- method,thatnonlinearinteractionsaresufficientforgen- sideredtheexpansionnecessarytotakeintoaccountfour erating power laws in wave spectra; moreover we show waveinteractionssothatweareconsistentwiththeorder that the Zakharov-Filonenko theory is completely con- of the “Zakharov equation”. sistent with the primitive equations of motion. We con- Adelicatepointinournumericalsimulationsisrelated sider a system of random waves localized in wave num- tothedissipationofenergyathighwavenumbers. Were- ber space and we show how nonlinearities “adjust” the markthatthisdissipationiscompletelyartificialsincewe spectruminagreementwiththeZakharovandFilonenko are dealing with a potential flow. Nevertheless we have prediction. Numerical work in the case of a forced and considered the dissipation phenomenon of the wave field dissipatedsystemhas been attempted by Willemsen [24] to be similar to the one that takes place in a turbulent using what sometimes are called the “Krasitskii equa- flow, i.e. that is mathematically expressed by a Lapla- tions” (see also [13]). In order to avoid the effects of cian that operates on the velocity. As is usually done external forcing, we considered the case of a freely de- in direct numerical simulations of box turbulent flows, cayingwavefield. Ifthesimulations,aswewillsee,show in order to increase the inertial range, we have used a the formation of a power law then the conjecture that higher order diffusive term. More explicitly on the right this power law is caused by geometrical features related handsideofequation(1)-(2),wehaveaddedrespectively to forcing and wave breaking must be excluded, since twoextraterms: −ν(−∇2)nψ and−µ(−∇2)mη,whereν forcingisabsentandwavebreakingcannotbetakeninto and µ represent an artificial viscosity coefficient and ∇2 account using the numerical method considered. From is the horizontal Laplacian. If n and m are greater than a physical point of view, the freely decaying case corre- 1 the viscosity is known as “hyperviscosity”. sponds to the evolution of a swell wave field. A generic It has to be noted that, at first sight, one would use a wavefield is consideredattime t=0 andit is allowedto veryhighpowerofthe Laplacianinordertoincreaseno- evolveinanaturalwayusingahighorderapproximation tablytheinertialrange,unfortunatelyveryhighvaluesof oftheEulerequations. Sincenumericalcomputationsare m and n could bring about the “bottleneck effect” [26], limited by the dimension of the grid considered, an arti- i.e. anaccumulationofenergyathighwavenumbersthat ficialdissipationisneededathighwavenumbersinorder coulddistortthe powerlawexpected[27]. Inournumer- to prevent accumulation of energy and a break down of icalsimulationsweusedν =µ=3×104 andn=m=8. the numerical code. The fluid is considered inviscid, ir- Thesevalueshavebeenselectedaftersometrialanderror rotational and incompressible. Under these conditions during the development of the numerical code: because the velocity potential φ(x,y,z,t) satisfies the Laplace’s of our limitated number of grid points, smaller values of equation everywhere in the fluid. The boundary condi- m and n, would obscure almost completely the inertial tions are such that the vertical velocity at the bottom is range. In our numerical simulations we did not impose zero and on the free surface the kinematic and dynamic any a dissipation at low wave numbers. boundary conditions are satisfied for the velocity poten- In order to prepare the initial wave field it is rea- tial ψ(x,y,t)=φ(x,y,η(x,y,t),t) (we assume that fluid sonable to consider a directional spectrum S(|k|,θ) = 2 P(|k|)G(θ). The directional spreading function G(θ) tions because, as time passes, energy is lost due to vis- used here is a cosine-squared function in which only the cosity, thus reducing the significant wave height of the firstlobe(relativetothedominantwavedirection)iscon- wave field and therefore the steepness. sidered: G(θ)= σ1cos2(cid:18)2πσθ(cid:19) if −σ ≤θ ≤σ (3) 101  0 else ) |  |k10-1 σ is a parameter that provides a measure of the direc- ( P tional spreading, i.e. as σ → 0, the waves become in- creasingly unidirectional. In our numerical simulations 10-3 t=0 h we selected the value of σ = π/2. We tried to avoid the t=0.1 h complete isotropic case in order to verify if the theory t=0.5 h still holds for intermediate values of the spreading. At 10-5 t=1 h the same time the selection of a large value of σ was motivated by the fact that recently it has been found 10-7 [28] that, for sufficiently narrow angle of spreading, the 0.01 0.1 |k| 1 Benjamin-Feir instability can be responsible for the formation of freak waves. As a consequence the nonlinear FIG.1. Wavespectra at different times. energy transfer could be slightly altered and some cor- rections to the prediction could be necessary (this very interesting topic is now under investigation and results willbereportedinadifferentpaper). P(|k|)ischosento EventhoughitisnotclearfromtheLog-Logrepresen- beanylocalizedspectrum. Wehaveperformednumerical tationinFig. 1,thereisadownshiftingofthepeakofthe simulations with a gaussianfunction or with a “chopped spectrumtowardslowerwavenumbers;asaconsequences JONSWAP”spectrum(aJONSWAPspectrumwitham- thesteepnesssubsequentlydecreasesovertime. Thetime plitudes equal to zero for frequencies greater than 1.5 scaleofthenonlinearenergytransferbecomeslargerand times the peak frequency) with random phases. For the larger. In Fig 2 we show the power spectrum of the sur- case of the gaussian function, wave numbers lower than faceelevationafter4hours(thesteepnesofthewavefield a selected threshold have been set to zero in order to is ε ≃ 0.07). In the same plot we show two power laws avoidextremely long largewaves. The velocitypotential ∼ k−2.5 and ∼ k−3: the first one seems to better fit the is then computed from the initial wave field using the data. In order to be completely sure that the numerical linear theory. Both gaussianand JONSWAP spectra led data are in agreement with the prediction of Zakharov to the same results in terms of the turbulent cascade. and Filonenko, we show in Fig. 3 compensated spectra Ourcomputationisperformedindimensionalunits;we with different compensation powers: z = 2.5 seems to haveselectedthe initialspectrumcenteredat0.1Hz,i.e. be the most plausible power law. Thus there seems to we are considering 10 seconds waves. The initial steep- be ample evidence from our numerical simulations that ness computed as ε = k0Hs/2 was chosen to be around the power law is in sufficiently good agreement with the 0.15 (Hs was computed as 4 times the standard devi- value predicted by the theory. ation of the wave field). The wave field was contained in a square grid (the resolution is 256×256) of length 103 L = 1417.6 meters. The time step considered was 1/50 the dominant frequency, i.e. ∆t=0.2 seconds. We have 102 performedournumericalsimulationsona400MhPC.In ) Fig. 1weshowtheevolutionofthewavepowerspectrum k|101 for different real times (t=0, 0.1, 0.5, 1 hours). We see (| P100 that,asexpected,thetailofthespectrumstartstogrow. This process seems to be quite fast: as is shown in the 10-1 figure after a few dominant wave periods some energy is already injected into high wave numbers. The process 10-2 k-2.5 of adjusting the power law to the “correct” one becomes k-3 then very slow, especially for low wave number. This 10-3 could be due to the frozen turbulent phenomenon [29], 10-4 i.e. a condition in which the energy fluxes towards high wave numbers are reduced because of the discretness of 0.01 0.1 |k| 1 the spectrum. Moreover decaying numerical simulations FIG. 2. Wave spectrum at t = 4 hours. A k−2.5 (dot- are very time consuming with respect to forced simula- ted-line) and a k−3 (dashed-line) power law are also plotted. 3 10-1 [3] V.E.ZakharovandN.Filonenko,SovietPhys.Dokl.,11, 881 (1967). z 10-2 [4] V.E.Zakharov,J.Appl.Mech.Tech.Phys.9,190(1968). | k [5] V.E. Zakharov, G.Falkovich and V. Lvov, Kolmogorov | |)10-3 spectra of turbulence I (Springer, Berlin, 1992). k [6] Y. Toba, J. Phys.Ocean., 29, 209 (1973). | P(10-4 z=2.5 [7] K. K. Kahama, J. Phys.Ocean. 11, 1503 (1981). [8] G. Z.Forristall, J. Geophys.Res. 86, 8075 (1981). z=3 [9] M.A.Donelan,J.HamiltonandW.H.Hui,Phil.Trans. 10-5 z=2 R. Soc. Lond. A 315, 509 (1985). [10] O.M. Phillips, J. Fluid Mech. 156, 505, (1985). 10-6 [11] D.ResioandW.Perrie, J.Fluid.Mech223,603(1991). [12] A. Pushkarev, D. Resio and V. Zakharov, submitted to 10-7 Phys. Rev.E (2001). 0.01 0.1 |k| 1 [13] V. E. Zakharov, Eur. J. Mech. B-Fluid 18 327, (1999). FIG. 3. Compensated wave spectra for different values [14] A. N. Pushkarev and V. E. Zakharov, Phys. Rev. Lett. of the compensation power: z = 2 (dahsed-line) z = 2.5 76, 3322, (1996). (solid-line) and z=3 (dotted-line). [15] DCai,A.J.Majda,D.W.McLaughlinandE.G.Tabak, Physica D 152, 551, (2001). [16] V. E. Zakharov, P. Guyenne, A. N. Pushkarev and F. After the pioneering work by Zakharovand Filonenko Dias, Physica D 152, 573, (2001). the kinetic wave theory has developed further, making [17] W.T.Tsai,D.K.P.Yue,Ann.Rev.FluidMech.28,249, available quantitative predictions for other physical ob- (1996). servablessuchasenergyfluxes,downshiftingofthepeak, [18] B. J. West, K. A. Brueckner and R. Janda, J. Geophys. energydissipationetc. Allthesequantitieswillbeexam- Res. 92, 11803, (1987). inedandresultswill be reportedinfuture papers. Other [19] D.G.DommermuthandD.K.P.Yue,J.FluidMech.184, questions naturally arise from our results: in HOS sim- 267 (1987). ulations the order of the computation depends on how [20] M. Tanaka, J. Fluid Mech. 442, 303 (2001). many terms are retained in the Taylor expansion; do [21] S. Yu. Annenkov and V.I. Shrira, J. Fluid Mech. 449, higher order terms influence the cascade process? Our 341 (2001). computation has been performed in a freely decaying [22] D.ClamondandJ.Grue,J.FluidMech.447,337(2001). case; could externalforcing (especially if anisotropic)in- [23] V. Zakharov, A. Dyachenko and O. Vasilyev, “New fluence the power law? And more, what would be the method for numerical simulation of a nonstationary po- influence of the water depth? These are all questions to tential flow of incompressible fluid with a free surface,” be answered in the near future. submitted for pubblication (2001). AcknowledgementsThisworkwassupportedbythe [24] J. Willemsen, Communication at the Theoretical devel- OfficeofNavalResearchofthe UnitedStates ofAmerica opments: two and three dimensional water waves, EU- (T. F. Swean, Jr.) and by the U.S. Army Engineer Re- Conference, Cambridge, UK,August (2001). search and Development Center. Consortium funds and [25] V.P. Krasitskii, J. Fluid Mech. 272, 1 (1994). Torino University funds (60 %) are also acknowledged. [26] G. Falkovich,Phys. of Fluids 6, 1411 (1994). M. O. was also supported by a Research Contract from [27] D. Biskamp, E. Schwarz, A. Celani, Phys.Rev. Lett. 81, 4855 (1998). the Università di Torino. [28] M.Onorato,A.R.OsborneandM.Serio,tobepublished in Phys. of Fluids, (2002). [29] A.N.Pushkarev and V.E.Zakharov, PhysicaD 135, 98 (2000). [30] L.E.Borgman,“TechniquesforComputerSimulationof OceanWaves”,inTopiccsinOceanPhysics, Proc. ofthe [1] A.Kolmogorov, C. R.Akad.Sci. SSSR30, 301 (1941). Int. School E. Fermi, Eds.A.R.OsborneandP.Malan- [2] O.M. Phillips, J. Fluid Mech 107, 426 (1958). otte Rizzoli, 373 (1980). 4