Turbulence of Capillary Waves on Shallow Water

Abstract

We consider the developed turbulence of capillary waves on shallow water. Analytic theory shows that an isotropic cascade spectrum is unstable with respect to small angular perturbations, in particular, to spontaneous breakdown of the reflection symmetry and generation of nonzero momentum. By computer modeling we show that indeed a random pumping, generating on average zero momentum, produces turbulence with a nonzero total momentum. A strongly anisotropic large-scale pumping produces turbulence whose degree of anisotropy decreases along a cascade. It tends to saturation in the inertial interval and then further decreases in the dissipation interval. Surprisingly, neither the direction of the total momentum nor the direction of the compensated spectrum anisotropy is locked by our square box preferred directions (side or diagonal) but fluctuate.

1. Introduction

Most of the ripples one observes on the puddles are capillary waves with wavelengths exceeding fluid depth, so they are one of the most ubiquitous waves in nature. Extensive literature on the turbulence of waves is summarized in two monographs [1,2]. As far as waves on the water surface are concerned, one finds a vast body of research on gravity waves and a substantial research on the capillary waves on a deep water. Yet surprisingly little is known about the turbulence of the capillary waves on a shallow water, except isotropic weak-turbulence spectrum and its small perturbations [1]. Needless to say that neither nature nor human activity generally provides us with nearly isotropic turbulence.

From a fundamental viewpoint it is of much interest to understand how anisotropy of forcing or geometry of the container impacts an anisotropy of developed turbulence at small scales. That problem has not been solved completely for capillary waves on a deep water either: Theory predicts that weak anisotropy imposed by forcing at large scales leads to more and more anisotropic turbulence at smaller and smaller scales, as the turbulence cascade develops [1,3]. On the other hand, initially anisotropic force-free turbulence was found numerically to undergo isotropization during turbulence decay [4], and also isotropization along the cascade was found for steady turbulence generated by a strongly anisotropic pumping [5].

This work is devoted to anisotropic turbulence of capillary waves on shallow water. The main research question that we pose is as follows: how anisotropy of the environment (forcing and/or container shape) influences the anisotropy of small-scale turbulence. Understanding far-from-equilibrium states of capillary waves on thin fluid layers in different geometries can be important, among other things, for an emerging field of liquid meta-materials—wave-driven matrices of vortices, akin to optical lattices [6]. Such wave-driven flows can give one an ability to control and separate active particles and chemicals in fluid layers in biological and engineering contexts, as well in controlled self-assembly [6,7]. The results we present below for capillary waves on shallow liquid shed some light on this phenomenon, which turns out quite generic.

We consider developed turbulence excited by a large-scale force producing zero wave momentum. We show that the energy cascade towards small scales spontaneously breaks the reflection symmetry $\mathbf{k}\to -\mathbf{k}$ and acquires non-zero momentum. The direction of this momentum fluctuates, yet even a long-time average in a fixed reference frame demonstrates substantial anisotropy. We also elucidate connection between anisotropy and interaction non-locality in k-space.

2. Methods

Here we present the model that describes small-amplitude capillary waves on shallow fluid and the analytical methods of its statistical analysis.

2.1. Basic Equations

Perturbations of the fluid surface can be described by the two scalar variables: the surface elevation $\eta (x,y,t)$ and the velocity potential $\mathsf{\Psi }(x,y,t)$. In the limit of small surface perturbations and shallow layer $\eta <h<\lambda$, the Hamiltonian has a very simple form [1]:

$$\begin{array}{c}\hfill \mathcal{H}=\frac{1}{2}\int \left({\sigma |\nabla \eta |}^2+\rho (h+\eta ){|\nabla \mathsf{\Psi }|}^2\right)dxdy\end{array}$$

(1)

Here h is the depth of an undisturbed fluid in cm, $\sigma$ is the surface tension in erg/cm ${}^2$ and $\rho$ is the density in g/cm ${}^3$. The two dynamical equations of motion are

$$\begin{array}{cc}& \frac{\partial \eta }{\partial t}=\frac{1}{\rho }\frac{\delta \mathcal{H}}{\delta \mathsf{\Psi }}=-\mathrm{div}(h+\eta )\nabla \mathsf{\Psi }=-h\Delta \mathsf{\Psi }-\mathrm{div}\eta \nabla \mathsf{\Psi },\hfill \end{array}$$

(2)

$$\begin{array}{cc}\hfill & \frac{\partial \mathsf{\Psi }}{\partial t}=-\frac{1}{\rho }\frac{\delta \mathcal{H}}{\delta \eta }=\frac{\sigma }{\rho }\Delta \eta -\frac{{|\nabla \mathsf{\Psi }|}^2}{2}.\hfill \end{array}$$

(3)

In the dimensionless variables, $\eta /h$, $\mathsf{\Psi }\sqrt{\rho /(\sigma h)}$, and $t\sqrt{\sigma /(\rho h^3)}$, we obtain:

$$\begin{array}{cc}& \mathcal{H}=\frac{1}{2}\int \left({|\nabla \eta |}^2+(1+\eta ){|\nabla \mathsf{\Psi }|}^2\right)dxdy\hfill \end{array}$$

(4)

$$\begin{array}{cc}\hfill & \frac{\partial \eta }{\partial t}=-\mathrm{div}\left((1+\eta )\nabla \mathsf{\Psi }\right),\hfill \end{array}$$

(5)

$$\begin{array}{cc}\hfill & \frac{\partial \mathsf{\Psi }}{\partial t}=-\frac{{|\nabla \mathsf{\Psi }|}^2}{2}+\Delta \eta .\hfill \end{array}$$

(6)

Wave amplitudes are defined as $a_{\mathbf{k}}=({\eta }_{\mathbf{k}}-i{\mathsf{\Psi }}_{\mathbf{k}})/\sqrt{2}$ and satisfy the equation of motion:

$$\begin{array}{c}\hfill \frac{\partial a_{\mathbf{k}}}{\partial t}=-ık^2{a}_{\mathbf{k}}-ı\int d\mathbf{q}{V}_{\mathbf{kqk}-\mathbf{q}}{a}_{\mathbf{q}}{a}_{\mathbf{k}-\mathbf{q}}-2ı\int d\mathbf{q}{V}_{\mathbf{qkq}-\mathbf{k}}^{*}{a}_{\mathbf{q}}{a}_{\mathbf{q}-\mathbf{k}}^{*},\end{array}$$

(7)

which corresponds to the Hamiltonian

$$\begin{array}{c}\hfill \mathcal{H}=\int {\omega }_k{|a_{\mathbf{k}}|}^{2}d\mathbf{k}+\frac{1}{2}\int \left(V_{\mathbf{qpk}}{a}_{\mathbf{q}}^{*}{a}_{\mathbf{p}}{a}_{\mathbf{k}}+\mathrm{c}.\mathrm{c}.\right)\delta (\mathbf{q}-\mathbf{p}-\mathbf{k})d\mathbf{q}d\mathbf{p}d\mathbf{k}\end{array}$$

(8)

where ${\omega }_k=k^2$ and $V_{\mathbf{qpk}}=p^2+k^2+(\mathbf{k}·\mathbf{p})$. In dimensional units ${\omega }_k=k^2\sqrt{\sigma h/\rho }$.

2.2. Simplest Properties

The system allows for 1D running wave of a permanent shape. Looking for $\eta (x-vt)$, $\mathsf{\Psi }(x-vt)$ and introducing $z=1+\eta$ one obtains $v{\mathsf{\Psi }}^{\prime }={(z{\mathsf{\Psi }}^{\prime })}^{\prime }$, which gives the mass flux constancy $z({\mathsf{\Psi }}^{\prime }-v)=V$, where the constant V is the solution parameter. Now we substitute that into $2z^{″}=-2v{\mathsf{\Psi }}^{\prime }+{\psi }^{\prime 2}=V^2{z}^{-2}-v^2$. Since we expect the profile inflection point where the surface is undisturbed, then $z^{″}=0$ when $z=1$ and $V=v$. Integrating once, we obtain

$$z^{\prime 2}/v^2=2C-z-1/z.$$

(9)

Integrating numerically this equation, one obtains the form of the solution for any v and $C\ge 1$. The wavelength is determined by the integral

$$\frac{\lambda }{2}=\frac{1}{v}\underset{C-\sqrt{C^2-1}}{\overset{C+\sqrt{C^2-1}}{\int }}\frac{d{z}^{\prime }}{\sqrt{2C-z^{\prime }-1/z^{\prime }}},$$

(10)

and must be commensurate with the box size, so that $L/\lambda =$ is an integer. For small amplitude waves, $z-1=\epsilon \ll 1$, which requires $C-1\ll 1$ we have the solution $\epsilon =\sqrt{C^2-1}sin(kx-k^{2}t)$, that is $v=k=2\pi /\lambda$. Existence of such solutions may hint that the system is integrable in one dimension, which can be expressed as a lack of thermalization, recurrence, etc. This is left for future studies.

3. Results

We now consider a set of many waves with small amplitudes. The only analytic results were obtained for the kinetic equation written for the occupation numbers defined according to $n_{\mathbf{k}}=\langle |a_{\mathbf{k}}{|}^2\rangle$. In the continuous limit the equation has the form [1,2]:

$$\begin{array}{cc}& \frac{\partial n_{\mathbf{k}}}{\partial t}=S{t}_{\mathbf{k}}=\int d{\mathbf{k}}_{1}d{\mathbf{k}}_2\left[U_{k12}(n_1{n}_{\mathbf{k}}-n_2{n}_{\mathbf{k}}-n_1{n}_2)-2U_{1k2}(n_2{n}_{\mathbf{k}}-n_1{n}_{\mathbf{k}}-n_1{n}_2)\right],\\ & U_{qpk}=\pi {|V_{\mathbf{qpk}}|}^2\delta ({\omega }_q-{\omega }_p-{\omega }_k)\delta (\mathbf{q}-\mathbf{p}-\mathbf{k}).\end{array}$$

(11)

It is valid when the interaction time $t_{NL}(k)\simeq n_{\mathbf{k}}/S{t}_{\mathbf{k}}$ exceeds the wave period, so that the dimensionless nonlinearity parameter ${ϵ}_k=1/{\omega }_k{t}_{NL}(k)={|V_{kkk}|}^2{n}_k{k}^d/{\omega }_k^2\ll 1$.

The kinetic Equation (11) conserves the energy $\int {\omega }_k{n}_{\mathbf{k}}d\mathbf{k}$ and the momentum $\mathbf{P}=\int \mathbf{k}{n}_{\mathbf{k}}d\mathbf{k}$. It has two known stationary isotropic solutions. The first one is the thermal equipartition ${\omega }_k{n}_k=$ const, so that in 2D the nonlinearity varies as a logarithm: ${ϵ}_k=ln(1/k)$. The second solution describes energy cascade in turbulence [1]:

$$n_k=k^{-4}.$$

(12)

or in dimensional variables

$$\langle |{\eta }_{\mathbf{k}}{|}^2\rangle =8{\epsilon }^{1/2}\sqrt{\pi }{(\rho h/\sigma )}^{1/4}{k}^{-4},$$

(13)

where $\epsilon$ is the energy flux. Such spectrum is expected to be realized in the inertial interval of wavenumbers between the pumping wavenumber $\kappa$ and the dissipation wavenumber K. For (12), ${ϵ}_k\propto k^{-2}$, so that if nonlinearity ${ϵ}_{\kappa }$ is small at the pumping scale, it is small everywhere. It is straightforward to check that the integrals in (11) converge on (12). Numerics demonstrate that the angle-averaged spectrum is indeed close to $k^{-4}$, as seen in Figure 1 below.

3.1. Analytic Theory of Angular Instability

To study the linear stability of the isotropic weak-turbulence spectrum (12), one takes $n_{\mathbf{k}}=k^{-4}+\delta n_{\mathbf{k}}$ and linearizes the kinetic equation with respect to $\delta n_{\mathbf{k}}$:

$$\begin{array}{cc}& \frac{\partial \delta n_{\mathbf{k}}}{\partial t}=\int d{\mathbf{k}}_{1}d{\mathbf{k}}_2\{U_{k12}\left[\delta n_1\left(k^{-4}-k_2^{-4}\right)+\delta n_2\left(k^{-4}-k_1^{-4}\right)-\delta n_{\mathbf{k}}\left(k_1^{-4}+k_2^{-4}\right)\right]\\ & -2U_{1k2}\left[\delta n_{\mathbf{k}}\left(k_1^{-4}-k_2^{-4}\right)+\delta n_2\left(k_1^{-4}-k^{-4}\right)-\delta n_1\left(k^{-4}+k_2^{-4}\right)\right]\}.\end{array}$$

(14)

Spectral analysis of this equation can be found in [8], here we employ the method of moments as the most compact way to establish the instability. We denote $\theta$ the polar angle of the vector $\mathbf{k}$ and define the moments of the relative perturbation $\delta n_{\mathbf{k}}/n_k=k^4\delta n_{\mathbf{k}}$:

$$M_m={\int }_{\kappa }^K{k}^{3}dk{\int }_0^{2\pi }\delta n_{\mathbf{k}}cos(m\theta )d\theta ,$$

(15)

Here $\kappa$ and K are respectively pumping and dissipation scales, setting the inertial interval of the turbulence cascade. We now substitute $V_{k12}$ and ${\omega }_k$ into (14) and integrate it to obtain the time derivatives of the moments. The most salient point is that convergence of the integrals are different for even and odd angular harmonics. For even m, two cancelations at small wavenumbers and one cancelation at large ones make the integrals convergent and thus independent of $\kappa \to 0$ and $K\to \infty$. Convergent integrals give $d{M}_m/dt$ proportional to m times a negative number (expressed via beta functions). The isotropic spectrum is thus linearly stable with respect to even harmonics. On the contrary, for odd harmonics $m=2j+1$, the cancelation proceed in powers of k rather than $k^2$, so that logarithmic divergencies remain. With logarithmic accuracy we obtain:

$$\frac{d{M}_m}{dt}={(-1)}^{j}mln(K/\kappa ).$$

(16)

We see that the isotropic spectrum must be linearly unstable with respect to double odd harmonics $m=1,5,\dots$

3.2. Numerical Modelling of Strongly Anisotropic Spectra

Analysis of substantially anisotropic turbulence is done by numerical solution of system () and () in double periodic domain. Partial derivatives are computed in Fourier space, using fast Fourier transform from FFTW library. Time advancement of $\mathsf{\Psi }$ and $\eta$, with right hand side as in () and (), is done with the 4th-order Runge-Kutta scheme. Next, to provide the small-scale dissipation, we set to zero the amplitudes of one third of highest modes; this common practice of dealiasing is done in Fourier space. At the same time and also in Fourier space, we force the lower 44 modes with $0<|k|<4$ by adding a noise, $f\xi k^{-2}\sqrt{\Delta t}$, to the Fourier images of $\mathsf{\Psi }$ and $\eta$, to provide forcing. Here f is the force amplitude and $\xi$ is the random noise white in time and uniformly distributed in the interval $[-1,1]$. The forcing is intended to model inhomogeneous pressure imposed on the surface by wind turbulence; such force must have zero total integral over space. Our forcing is substantially anisotropic, yet statistically invariant with respect to the transformation $\mathbf{k}\to -\mathbf{k}$, so that on average it does not produce any momentum. On the other hand, we have shown above that unstable harmonics include $m=1$ which has non-zero momentum. The time steps are chosen based on spatial resolution, $\Delta t=0.2{(\Delta x)}^2$, to resolve the frequencies of the highest modes, ${\omega }_k=k^2$. This is an analog of CFL condition for our dispersion relation, with empirical coefficient 0.2 determined in test simulations.

Figure 1 shows that the angle-averaged spectrum is indeed close to $k^{-4}$, as predicted by (12,13). We see that weak turbulence work well even for a sufficiently strong forcing: $f=0.005$ corresponds to ${\eta }_{rms}\simeq 0.2$.

The most interesting questions are then as follows:

• Is there a spontaneous breaking of symmetry $\mathbf{k}\to -\mathbf{k}$ and generation of nonzero total momentum of the wave system?
• How does the degree of anisotropy behave as one goes along the energy the cascade to higher k?

Figure 2 and Figure 3 answer the first question. Figure 2 presents a long run which includes the spectrum formation. We see that nonzero momentum spontaneously appears and persists. Its direction wanders around, but here is mostly oriented along the box sides ($\theta /\pi \simeq 0.5$) or diagonal ($\theta /\pi \simeq 0.75,0.25$).

Figure 3 presents the results of 10 runs starting with the same initial data of established spectrum with randomized phases. The runs are subjected to different random force realizations. For turbulence, the momentum ${\sum }_{\mathbf{k}}\mathbf{k}{n}_{\mathbf{k}}(t)$ is mostly determined by the lowest wavenumbers. To characterize the overall degree of turbulence anisotropy, we define the vector ${\sum }_{\mathbf{k}}\mathbf{k}{k}^4{n}_{\mathbf{k}}$. We find that its modulus is of order unity in all cases, that is all the spectra are substantially anisotropic.

Figure 4, Figure 5 and Figure 6 answer the second question (on the k-dependence of the degree of anisotropy). Here we see that pumping-connected strongly anisotropic spectrum at low k undergoes partial isotropization as one proceeds along the cascade and tend to saturate at moderate anisotropy in the inertial interval. Figure 4 shows the average over a long run (about 25,000 dimensionless times) presented in Figure 2; we have checked that averaging in the fixed reference frame gives about the same spectra as averaged in the local frame directed along the instantaneous anisotropy vector ${\sum }_{\mathbf{k}}{k}^4{n}_{\mathbf{k}}(t)$.

One obtains similar results averaging over in the local frame over 10 relatively short runs (4000 dimensionless time), all started from the same state with randomized phases, as shown in Figure 5, which is counterpart to Figure 3.

Note how the angular form of the spectrum in the inertial interval approaches the first harmonic in agreement with the theoretical prediction (16). Further, we compute the angular harmonics of the spectra

$$F_m(k)={\int }_0^{2\pi }{n}_{\mathbf{k}}cos(m\theta )d\theta ,$$

(17)

Figure 6 shows the normalized first two angular harmonics for short and long averages. The data are too noisy for harmonics with $m>2$ to be reliable. One sees that the second harmonic decays with k much faster than the first one. Comparison of the dotted and dashed lines for the first harmonic shows that anisotropy in a static reference frame tend to be smaller, apparently because the direction of the anisotropy vector fluctuates in time. Yet the static average depends on k only weakly, which shows that anisotropy persist in the inertial interval even on such a long timescale (over 25,000 dimensionless times).

For interpretation of Figure 6, it is instructive to recall the difference in k-dependence of angular perturbations to thermal equilibrium and to turbulence. Wave system with a nonzero total momentum (drift) has thermal equilibrium $n_{\mathbf{k}}\propto {[{\omega }_k-(\mathbf{k}·\mathbf{v})]}^{-1}$, where $\mathbf{v}$ is the drift velocity. Therefore, the angular harmonics for weakly anisotropic spectra behave as follows:

$$\frac{F_m(k)}{F_0(k)}\propto {\left(\frac{kv}{{\omega }_k}\right)}^m\propto k^{-m}.$$

(18)

On the contrary, stationary anisotropic corrections to turbulence spectra cannot be obtained from the isotropic spectrum by the Galilean transform ${\omega }_k\to {\omega }_k-(\mathbf{k}·\mathbf{v})$ [1,9]. The anisotropic correction carrying small momentum flux $\mathbf{R}$ from pumping to damping must be proportional to the dimensionless ratio of the fluxes: $F_m(k)/F_0(k)\propto (\mathbf{R}·\mathbf{k}){\omega }_k/\epsilon k^2\propto k$, which grows with k. Since our low-k pumping does not produce any momentum, we do not expect any spectral momentum flux and such behavior. Indeed, Figure 6 does not show any growth with k in relative amplitude of harmonics. Surprisingly, the short-run averages behave according to (18) in the inertial interval (red and blue solid lines in Figure 6) despite this being turbulence rather than thermal equilibrium. On the other hand, long-run averages of the first harmonic do not follow these laws but rather tend to saturate (broken and dotted red lines in Figure 6.

4. Discussion

We first briefly discuss the applicability limits of our specific findings and then list them and present the possible wider implications. Specifically, the equations that we solve analytically and numerically describe waves with the wavelengths longer than the fluid depth yet shorter than $\sqrt{\sigma /\rho g}$, which for most fluids ranges between one and several centimeters. We described developed turbulence, that is assume a relatively wide transparency window, that is the interval of wavelength where the effects of viscous dissipation can be neglected relative to the nonlinear interaction between waves. For clean fluids (with viscosities of order $0.1÷0.01$ cm ${}^2$/sec) even for weak nonlinearity $\eta /h<0.1$, dissipation effects start to be important only for sum-millimeter wavelengths. However, special care is needed in applications to particle-laden surfaces, where both the surface tension and the dissipation rate may vary significantly.

Classical weak-turbulence theory gives an isotropic spectrum with local interaction and local cascade in k-space. However, the isotropic spectrum is unstable with respect to the anisotropic angular perturbations of double odd orders, $m=4n+1$, which interact non-locally both with smallest and largest wavenumbers. That instability must lead, in particular, to the spontaneous breakdown of the reflection symmetry and generation of a nonzero momentum of the wave system.

Direct numerical simulations support the prediction and show substantially anisotropic spectra without reflection symmetry. We apply an anisotropic large-scale random pumping having the square box symmetry, thus producing on average zero momentum. In all runs, the spectra spontaneously acquire nonzero momentum. Directions of the momentum and anisotropy vectors wander with time. It may well be that averages on astronomically large timescales produce isotropic spectra; we have not been able to see any tendency towards that even on the largest timescales of tens of thousands dimensionless times.

The degree of spectrum anisotropy decreases along the cascade and tend to saturate in the inertial interval, so that the dimensionless parameter of anisotropy is of order unity, that is the anisotropy is neither weak nor strong. That may have a wider implications. For example, as mentioned in the Introduction, isotropization of small scales for decaying turbulence [4] and for a steady cascade [5] was found for capillary waves on deep water, despite theoretical prediction for linear instability with respect to angular perturbations [1,3]. It is likely that this seeming contradiction can be explained similarly to what we observe: stationary spectrum deep in the inertial has a moderate degree of anisotropy.