2.1 Experimental set-up

In the present work, we use the experimental results presented in Scolan et al. (Reference Scolan, Ermanyuk and Dauxois2013) and also additional experimental data obtained with the experimental set-up sketched in figure 1(a). A Cartesian coordinate system is introduced, with the horizontal $x$ and vertical $z$ axes located in the vertical midplane of the test tank. The $y$ axis is perpendicular to the $(x,z)$ plane and directed from the observer/camera towards the tank. The rectangular test tank of size $800\times 170\times 425$  mm $^{3}$ is filled with a salt-stratified fluid using the conventional double-bucket technique. The density distribution as a function of the vertical coordinate $z$ is measured prior and after experiments by a conductivity probe driven by a vertical traverse mechanism. The value of the buoyancy frequency $N=[(-g/\bar{{\it\rho}})(\text{d}{\it\rho}/\text{d}z)]^{1/2}$ is inferred from the linear fit to the measured density profile. An example of such a profile is shown on figure 1(b). A sliding sloping wall, inclined at the angle ${\it\alpha}$ , is carefully inserted into the fluid once the filling procedure is over. The sloping wall delimits a trapezoidal fluid domain of length  $L$ (measured along the bottom) and depth $H$ .

The input forcing is introduced into the system by an internal wave generator (Gostiaux et al. Reference Gostiaux, Didelle, Mercier and Dauxois2007; Mercier et al. Reference Mercier, Martinand, Mathur, Gostiaux, Peacock and Dauxois2010; Joubaud et al. Reference Joubaud, Munroe, Odier and Dauxois2012; Scolan et al. Reference Scolan, Ermanyuk and Dauxois2013) with the time-dependent vertical profile given by

where $a$ and ${\it\omega}_{0}$ are the amplitude and frequency of oscillations, respectively. The profile is reproduced in discrete form by the motion of a stack of $47$ horizontal plates driven by the rotation of a vertical shaft. The stack of moving plates is confined between lateral walls of a vertical box-like support frame. Each lateral wall of the frame is $15$  mm thick. Therefore, the inner size of the frame, $W^{\prime }$ , which coincides with the width of the moving plates, is $30$  mm smaller than the width of the tank $W$ . Thus, the forcing is applied to the stratified fluid with a width $W^{\prime }$ representing $82\,\%$ of the whole width $W$ .

Fluid motion is recorded by a computer-controlled video AVT (Allied Vision Technologies) Stingray camera with CCD matrix of $1388\times 1038$ pixels or Pike camera with CCD matrix of $2452\times 2054$ . The camera is placed at a distance between 1750– $2300$  mm from the test tank. Given the relatively small size of the working section, no correction for the parallax errors is applied. Two methods are used to measure the parameters of the fluid motion. The perturbations of the density gradient are measured with synthetic schlieren technique (Sutherland et al. Reference Sutherland, Dalziel, Hughes and Linden1999; Dalziel, Hughes & Sutherland Reference Dalziel, Hughes and Sutherland2000) from apparent distortions of a background random dot pattern observed through the stratified fluid. In some cases, standard PIV technique has alternatively been used for velocity measurements. The fluid is seeded with light-reflecting particles of typical size $8~{\rm\mu}\text{m}$ and density $1.1~\text{kg}~\text{m}^{-3}$ . The sedimentation velocity of particles is found to be sufficiently low. Its effect on the results of velocity measurements is negligibly small. The test section is illuminated by a vertical laser sheet coming through the transparent side or bottom of the test tank. The spanwise position of the sheet could be varied to assess 3D effects.

For the synthetic schlieren technique, the camera takes two frames per second. For the PIV measurements, the camera operates in two-frame burst mode with an interval between the pair of images within a burst set to $0.125$  s and a global burst cycle period of $0.5$  s. Conversion of images into velocity fields is performed using a cross-correlation PIV algorithm with subpixel resolution (Fincham & Delerce Reference Fincham and Delerce2000). Thus, with two velocity fields per second, i.e. typically around $20$ fields per wave period, we can resolve the significant frequency components of the signal. In addition, the spatial resolution for the velocity field is about 3 mm in each direction and is found sufficient to resolve fine details of the wave field.

2.2 Numerical set-up

The numerical simulation of wave attractors faces two major challenges. First, the fluid motion is highly nonlinear, and accounting for nonlinear interactions is crucial for the dynamics, even for weak interactions. Second, the Schmidt number ( $Sc$ ), defined as the ratio of water kinematic viscosity and salt diffusivity, is close to 700 in a salt-stratified fluid. So the spatial scale of the density perturbations can be much smaller than the scale of the velocity perturbations and more demanding in terms of spatial resolution. The scale under which no scalar gradient remains because of diffusion effect is very small and smaller than the usual Kolmogorov scale $L_{K}$ used as a mesh criterion for direct simulations where $L_{K}=Re^{-3/4}\times L$ with $L$ the integral scale and $Re$ the large-scale Reynolds number. Strictly speaking, this small scale where scalar diffusion takes place can be estimated with the Batchelor scale $L_{B}=L_{K}\times Sc^{-1/2}$ (Batchelor Reference Batchelor1959; Buch & Dahm Reference Buch and Dahm1996; Rahmani, Seymour & Lawrence Reference Rahmani, Seymour and Lawrence2014).

The first numerical simulation of the formation of internal wave attractors (Grisouard et al. Reference Grisouard, Staquet and Pairaud2008) used finite volume method, one of the most popular state-of-the-art tools in computational fluid dynamics. However, as noted by these authors, the numerical simulation could only reproduce the dynamics of the attractors for Schmidt number less than 100. In addition, the discretization of convective terms produces a numerical viscosity, which acts similarly to a real viscosity. This numerical artefact blurs the fine-scale structures arising due to the high value of the Schmidt number and it may also introduce substantial errors in calculations of the long-term dynamics of attractors (at time scales of tens or hundreds of periods). However, for the current study, the long time intervals are precisely of particular interest to study the growth, saturation and cascades of secondary instabilities.

As an alternative, the Galerkin method represents a highly efficient approach to tackle the nonlinear interactions if a full system of basis functions subject to boundary conditions can be used. A pseudospectral approach offers an efficient realization of the Galerkin method. In Jouve & Ogilvie (Reference Jouve and Ogilvie2014), this approach has been used to study nonlinear interactions in attractors of inertial waves. In the present work we have chosen the method of spectral elements which combines the accuracy and high resolution of spectral methods with geometric flexibility of finite element methods, and what is particularly suitable for simulation of long-term evolution of fine-scale flows in globally forced geophysical systems (Favier et al. Reference Favier, Grannan, Le Bars and Aurnou2015). The computational domain is divided into a finite number of quadrilateral (in two dimensions) or hexahedral (in three dimensions) elements and a Galerkin approach is applied to each element. For the numerical realization, we use the open code Nek5000, developed by Paul Fischer and colleagues (Fischer & Ronquist Reference Fischer and Ronquist1994; Fischer Reference Fischer1997; Fischer & Mullen Reference Fischer and Mullen2001). In each element, the Lagrange polynomial decomposition is used and applied at Gauss–Lobatto–Legendre points for the sake of stability (to avoid ill-conditioning). The full resulting mesh, consisting of the elements and the Gauss–Lobatto–Legendre points, is highly non-uniform, adding a bit more complexity in post-processing data treatment. However, the efficiency of the code fully justifies such a nuisance. In other words, the approach used weighted residual techniques employing tensor-product polynomial bases. In addition other benefits, it allows ‘analytical’ computation of the derivatives through matrix–matrix products or matrix–matrix-based evaluation

in which $u$ represents any of the unknown variables, for instance, velocity components, density, etc., and ${\it\psi}$ is a Lagrangian interpolant through the Gauss–Lobatto–Legendre points.

The implementation of a high-order filter in Nek5000 allows us to stabilize the method for convection-dominated flows (Fischer & Mullen Reference Fischer and Mullen2001). The time advancement is based on second-order semi-implicit operator splitting methods and stable backward-difference scheme. The additive overlapping Schwarz method is used as a pre-conditioner (Fischer Reference Fischer1997).

The most difficult regions from a computational viewpoint are located in the vicinity of the rigid walls since, in these regions, intense folding of high-gradient density layers may occur, especially in the case of a high-amplitude wave motion. For this reason, a non-uniform element mesh in the near-wall regions is preferable for simulations of the nonlinear dynamics of wave attractors. Typically, we have used meshes with up to 0.5 million elements, with eighth- to tenth-order polynomial decomposition within each element.

The full system of equations being solved consists of the Navier–Stokes equation in the Boussinesq approximation

the continuity equation

and the equation for the transport of salt

where ${\it\rho}_{m}$ is the density of solution with the constant minimal reference salinity, ${\it\rho}_{s}$ the density perturbation in a unit volume due to local salinity (the full density is ${\it\rho}={\it\rho}_{m}+{\it\rho}_{s}$ ), ${\it\mu}$ the dynamical viscosity and ${\it\lambda}$ the diffusivity of salt. Dynamical viscosity and diffusivity are assumed to be constant. We impose no-slip boundary conditions on the rigid surfaces and stress-free condition on the upper surface. The boundary conditions on ${\it\rho}_{s}$ are isolation: $\partial {\it\rho}_{s}/\partial n=0$ where $n$ is the normal to the wall. Forcing is applied at the vertical wall by prescribing the profile of the horizontal velocity which reproduces the motion of the generator (2.1) in the vertical direction and takes into account the difference between the width of the generator $W^{\prime }$ and the width of the tank $W$ . The transverse profile of the forcing is prescribed by stepwise or piecewise linear function. At the central segment of the width $W^{\prime }$ , the horizontal velocity is uniformly distributed in the transverse direction in both cases. At the side segments of width $(W-W^{\prime })/2$ , we prescribe either (i) zero velocity or (ii) a linear decrease of the velocity from the uniform value to zero. These two versions of transverse profiles of forcing give essentially the same results, with a small quantitative difference in the second version where forcing gives a higher horizontal impulse to the system. Both versions can be implemented in calculations, but the second version of the forcing is found to be much less expensive computationally and the major part of computations was performed with it. The efficiency of the generator used in experiments is not 100 %. For this reason, the amplitude of forcing in computations with the piecewise linear case is intentionally reduced by 10 % compared to the experimental value, thereby providing a good match between numerical and experimental results.

The comparison between numerical and experimental results is typically complicated by the presence of a thin mixed layer close to the free surface of the stratified fluid in the test tank. Because of this layer, the internal wave beam is not reflected precisely at the free surface. Instead the wave beam undergoes a complex reflection, partially at the interface between the mixed layer and the linearly stratified fluid and partially at the free surface. This complex reflection affects the shape of the attractor (Grisouard et al. Reference Grisouard, Staquet and Pairaud2008) and the shape of envelopes of wave motions in the wave beams. To take this layer into account in the numerical model, we introduce a model density distribution where the full depth is denoted $H$ as in experiments, the thickness of the thin mixed upper layer is ${\it\delta}$ , and the depth of the linearly stratified fluid is $H^{\prime }=H-{\it\delta}$ .

3.1 A stable attractor

Figure 2.

Experimental (a) and numerical (b) snapshots of the horizontal density gradient at $t=50T_{0}$ . The amplitude of the wave maker is $a=2~\text{mm}$ for the experiment and $a=1.8~\text{mm}$ for the simulation. Both attractors are stable. Note that the shade (colour online) scale is the same in both panels. The small black quadrilateral defines the acquisition region used for computing the time–frequency spectrum presented in figure 4. The wave frequency is ${\it\omega}_{0}/N=0.62\pm 0.01$ . Experimental (c) and numerical (d) horizontal density gradients as a function of the time, for a point located on the most energetic branch within the black trapezoid depicted in panels (a, b). In the calculation we take a piecewise linear approximation of the experimental density profile, with the lower layer of depth $H^{\prime }=308~\text{mm}$ and buoyancy frequency $N$ , and the upper layer of depth ${\it\delta}=18~\text{mm}$ with a density gradient 8 times smaller. The total depth of fluid $H=H^{\prime }+{\it\delta}=326~\text{mm}$ .

As already mentioned in the introduction, the dynamics of stable attractors has been studied in great detail before, both in stratified and rotating fluids. It is well established that the typical thickness of the attractor scales as $1/3$ power of the Ekman/inverse Stokes number in rotating/stratified fluid (Rieutord et al. Reference Rieutord, Georgeot and Valdettaro2001; Ogilvie Reference Ogilvie2005; Grisouard et al. Reference Grisouard, Staquet and Pairaud2008; Hazewinkel et al. Reference Hazewinkel, van Breevoort, Dalziel and Maas2008). The evolution of the spatial wave spectrum during the transient regime after the start-up of the oscillations toward an equilibrium shape is also well-known (Grisouard et al. Reference Grisouard, Staquet and Pairaud2008; Hazewinkel et al. Reference Hazewinkel, van Breevoort, Dalziel and Maas2008). This dynamics is well-reproduced in our calculations be they 2D or 3D. In addition, we observe a very good quantitative agreement in our 3D simulations once we introduce a small correction for non-perfect efficiency of the wave generator.

Typical snapshots of the computed and measured fields of the horizontal density gradient ${\it\rho}_{x}$ are shown in figure 2, emphasizing a very good qualitative and quantitative agreement between the numerical and experimental data. This good agreement is further confirmed with the time series of the horizontal density gradient at a point located in the most energetic branch of the attractor as illustrated in figure 2. Note that the numerical value of the density gradient is computed at the vertical midplane $xOz$ of the test tank. With this comparison, we tacitly assume that the flow is approximately 2D and the standard schlieren technique is applicable (Dalziel et al. Reference Dalziel, Hughes and Sutherland2000). The effects of three-dimensionality are discussed further in the next section.

3.2 An unstable attractor

As the amplitude of oscillations of the wave generator increases, the attractor becomes unstable. In what follows, we describe the comparison of the experimental and numerical results for the onset of the TRI. Two-dimensional simulations were also performed and were found to predict the instability to occur at the forcing amplitudes roughly two times smaller than the experimental ones. This is due to the absence of lateral walls in the 2D simulations. These walls introduce a significative dissipation in 3D simulations, as demonstrated in § 4.5. In view of this quantitative discrepancy between 2D and 3D simulations, in what follows we use the results of 3D calculations for comparison with the experiments. The experimental and numerical snapshots of the horizontal density gradient field are presented in figure 3. The development of TRI is clearly seen in the most energetic branch of the attractor. A very good quantitative and qualitative agreement is again observed between the experimental and numerical wave fields. The experimental and numerical time series at a point in the most energetic branch of the attractor show also an excellent correspondence.

Figure 3.

Experimental (a) and numerical (b) snapshots of the horizontal density gradient at $t=50T_{0}$ . The amplitude of the wave maker is $a=2.5~\text{mm}$ for the experiment and $a=2.4~\text{mm}$ for the numerical simulation. The wave frequency is ${\it\omega}_{0}/N=0.62\pm 0.01$ . Both attractors are unstable and the instability appears first on the most energetic branch. The small black quadrilateral defines the acquisition region used for computing the time–frequency spectrum presented in figure 4. Experimental (c) and numerical (d) horizontal density gradients as a function of the time, for a point located on the most energetic branch within the black trapezoid depicted in panels (a,b). The density stratification is the same as described in the caption to figure 2.

Let us now focus on the development of the instability in more detail to see whether computations reproduce the experimentally observed triads in temporal and spatial domains. We first consider the evolution of the instability in the temporal domain. The development of the frequency spectrum of wave motion over time is presented in figure 4. The time–frequency diagrams are calculated from numerical and experimental data for points as in Bourget et al. (Reference Bourget, Dauxois, Joubaud and Odier2013), with the formula

where $h$ is an Hamming window and $r$ the component of the velocity field along $x$ or $z$ , i.e.  $S_{u}$ or $S_{w}$ . The calculations are performed with the Matlab toolbox described in Flandrin (Reference Flandrin1999). The appropriate choice of the length of the Hamming window allows to tune the resolution in space and time. The experimental and numerical spectra are averaged over the analysing area shown in figure 3. It can be seen that the numerical and experimental spectra agree qualitatively and quantitatively. The signal is initially entirely dominated by the forcing frequency ${\it\Omega}_{0}=0.62$ (in which ${\it\Omega}_{0}={\it\omega}_{0}/N$ is the non-dimensional frequency of oscillation): it corresponds precisely to the primary (carrier) wave. Then oscillations with frequencies ${\it\Omega}_{1}$ and ${\it\Omega}_{2}$ , which correspond to two secondary waves generated by TRI, slowly develop with time. At $t=50T_{0}$ , these frequencies are ${\it\Omega}_{1}=0.24$ and ${\it\Omega}_{2}=0.38$ . They satisfy the frequency conditions for the triadic resonance

In addition, one can see also two peaks ${\it\Omega}_{3}=0.86$ and ${\it\Omega}_{4}=1.00$ satisfying differential conditions

To verify the fulfillment of the condition for triadic resonance in space, we apply the Hilbert transform technique (Mercier, Garnier & Dauxois Reference Mercier, Garnier and Dauxois2008; Bourget et al. Reference Bourget, Dauxois, Joubaud and Odier2013) to the numerically simulated data in the same way as is done in experiments (Bourget et al. Reference Bourget, Dauxois, Joubaud and Odier2013): the experimental signal is demodulated with the Fourier transform, filtered around the three frequencies of interest, ${\it\Omega}_{0}$ , ${\it\Omega}_{1}$ and ${\it\Omega}_{2}$ , and reconstructed back in real space using the inverse Fourier transform. As shown in figure 5 of Scolan et al. (Reference Scolan, Ermanyuk and Dauxois2013), it is possible to compute amplitude and phase of each component. The latter appear as patterns of stripes, corresponding to a fixed moment of time. The wave vectors can be derived by differentiating these phases along the $x$ and $z$ directions.

Figure 4.

Experimental (a) and numerical (b) time–frequency spectra of data presented in figure 3. Both have been computed with the same signal processing parameters, over the same area, located on the most energetic branch and depicted with a black trapezoid in figure 3.

To quantify the wave vectors involved into the triadic resonance in the numerical data, we construct probability density functions (PDFs) for the components of wave vectors. Using these PDFs for the components of the wave vector corresponding to the primary wave oscillating at ${\it\Omega}_{0}$ , we can estimate $\boldsymbol{k}_{0}=\{-63.5,-80.1\}~\text{m}^{-1}$ . The same procedure for the secondary wave oscillating at ${\it\Omega}_{2}$ gives $\boldsymbol{k}_{2}=\{-104.8,-242.5\}~\text{m}^{-1}$ . The phase pattern is slightly more complicated for the component oscillating at ${\it\Omega}_{1}$ as can be noticed in figure 5 of Scolan et al. (Reference Scolan, Ermanyuk and Dauxois2013) where lines of equal phase for this frequency are not completely straight and there is a larger error bar on the measurement of the wave vector. Numerically, it is possible to investigate further in the focusing branch of the attractor thanks to the PDF analysis: one realizes that there are two families of collinear wave vectors, which yield multi-peaked PDF of the wave vector components, where the main peak corresponds to $\boldsymbol{k}_{1}=\{46.5,174.3\}~\text{m}^{-1}$ . These values of the wave vector components closely match the spatial condition of the triadic resonance. The physical interpretation of the other peaks present in the PDF of the wave vector components for the wave field filtered at ${\it\Omega}_{1}$ is unclear.

Figure 5.

Experimental (a) and numerical (b) triads. The three wave vectors are measured from the phase of the waves after Hilbert filtering around ${\it\omega}_{0}$ , ${\it\omega}_{1}$ and ${\it\omega}_{2}$ . The error bars are given by of ellipses at the arrow tips.

The triangles of the wave vectors obtained in experiments and simulations are shown in figure 5. Again, we observe a good quantitative agreement between the experimental and numerical results.

4.2 Two-dimensional visualizations in vertical plane $xz$ at different transversal locations

In Scolan et al. (Reference Scolan, Ermanyuk and Dauxois2013), the experiments presented have been performed using the synthetic schlieren technique which assumes that the flow is bidimensional. Indeed, in this technique, the image results from the integration of the field along the transverse direction and the ray of light is expected to cross the tank through always the same density perturbation. This assumption seems to be reasonable because the experimental set-up has been designed on purpose to ensure the two-dimensionality of the flow. Indeed, the tank is narrow in the transverse $y$ direction and forcing is applied via the wave maker through almost all the width of the tank so that the flow is expected to be 2D, except very close to both vertical lateral walls where the viscosity plays an important role.

To check this assumption, we use experimental standard PIV technique for several vertical sections at different constant $y$ positions of the light sheet and thus compare flow properties between the different sections. Experiments have been performed with $L=435$  mm, $H=300$  mm and ${\it\alpha}=25^{\circ }$ , a geometry close to experiments performed in Scolan et al. (Reference Scolan, Ermanyuk and Dauxois2013). Different vertical sections in the transverse direction have been illuminated with a vertical laser sheet, coming from the bottom of the tank, which is transparent. A $45^{\circ }$ inclined mirror, located below and as long as the tank, transforms a horizontal laser sheet into a vertical one. By simply translating the mirror in the transverse direction, one can illuminate the different test sections. Three series of experiments have been performed, each one with a different amplitude of the wave maker: $a=1.5$ , $3$ and $5$  mm. For each series, between 9 and 12 vertical sections have been illuminated by simply translating the mirror in the transverse direction.

Measurements have been carried out in the stationary state of the attractor. For the lowest amplitude experiment, there is no instability of the attractor: the wave maker was stopped only at the end of the full series, once all sections have been illuminated. For the two other series, the instability takes place after a long transient: consequently, the wave maker has been stopped after performing measurements at each transversal location to avoid mixing effects that would modify the stratification and thus the attractor itself. For all series, measurements in the central section of the tank were performed twice, once at the beginning, once at the end. This gives an estimate of the errors made on these measurements and also gives the possibility to check that, after a whole series, the attractor is really unchanged. To compare the different sections of the same series, we focused only on the most energetic branch of the attractor (the one connecting the slope to the free surface). Transversal profiles are extracted through this branch, always at the same location. As the horizontal and vertical velocities are measured, it is possible to calculate the velocity along this branch inclined at the angle ${\it\theta}=\arcsin {\it\Omega}_{0}$ using the formula $v_{{\it\xi}}=-v_{x}\cos {\it\theta}+v_{z}\sin {\it\theta}$ . The maximum value corresponds to the centre of the branch, ${\it\eta}=0$ , where ${\it\eta}$ is the coordinate transversal to the branch.

Numerically, with 3D simulation, it is possible to get velocity profiles in the transverse $y$ direction. Two 3D simulations have been performed, with exactly the same geometrical parameters as in Scolan et al. (Reference Scolan, Ermanyuk and Dauxois2013), and with two different amplitudes of the wave maker: $a=1.6$ and $2.2$  mm. Velocity profiles in the transverse $y$ direction were extracted from the most energetic branch and the velocity $v_{{\it\xi}}$ along this branch computed as for experiments.

It is important to note that the transverse resolution is much higher in the numerical data ( $512$ points for $170$  mm) than in the PIV experimental ones (between $9$ and $12$ points for the same $170$  mm). Moreover, the different points obtained experimentally come from different individual experiments performed at the same input parameters but not from the same unique experiment. Thus, it is difficult to synchronize the different points and we are limited to look at amplitudes of the filtered signals at ${\it\omega}_{0}$ , ${\it\omega}_{1}$ and ${\it\omega}_{2}$ . In contrast, various quantities can be examined numerically such as the raw signal, the whole signal filtered around ${\it\omega}_{0}$ , ${\it\omega}_{1}$ and ${\it\omega}_{2}$ and the mean flow generated by the attractor. Different videos, available in the supplementary information, show these different fields.

4.3 Velocity profiles in transverse $y$ direction

Figure 7.

(a) Raw velocity profiles as a function of the transversal direction $y$ for the numerical simulation with $a=1.6~\text{mm}$ . Thick solid line and thin solid line: maximum ( $\max v_{{\it\xi}}$ ) and opposite of the minimum ( $-\!\min v_{{\it\xi}}$ ) of the velocity $v_{{\it\xi}}$ for each $y$ over the time history. Dashed line: maximum of the absolute value of the transversal velocity $v_{y}$ over the time history for each  $y$ . (b) Modulus of $v_{{\it\xi}}$ filtered over $20T_{0}$ around ${\it\omega}_{0}$ (thick dashed line), for the numerical simulation with $a=1.6~\text{mm}$ . Modulus of $v_{{\it\xi}}$ filtered over $20T_{0}$ around ${\it\omega}_{0}$ (thick solid line), ${\it\omega}_{1}$ (thin solid line) and ${\it\omega}_{2}$ (thin dashed line), for the numerical simulation with $a=2.2~\text{mm}$ . (c) Normalized profiles of $v_{{\it\xi}}$ filtered around ${\it\omega}_{0}$ for the numerical simulation with $a=1.6~\text{mm}$ (thick dashed blue line) and $a=2.2~\text{mm}$ (thin solid red line). Points were obtained using PIV experiments: $a=1.5~\text{mm}$ (blue diamond), $a=3~\text{mm}$ (red square) and $a=5~\text{mm}$ (black circle).

Figure 7(a) presents the different velocity profiles for the numerical simulation with $a=1.6~\text{mm}$ . For this amplitude, the attractor is stable. The raw velocity profiles for $v_{{\it\xi}}$ are, as expected, very close to the bidimensionality, except into the boundary layers where the velocity goes to zero. Moreover, it appears clearly that the amplitude of the transverse velocity $v_{y}$ (represented by a dashed line) is much smaller than the velocity $v_{{\it\xi}}$ along the branch. This is a clear confirmation that the attractor generated by the wave maker is 2D to a good approximation. The difference between the raw velocity magnitudes measured in positive and negative ${\it\xi}$ directions comes from the presence of slow near-wall secondary currents superposed with the fundamental monochromatic wave motion in the wave beams as discussed in more detail in § 4.4.

To extract only the component of the wave field oscillating at the fundamental frequency, one can perform Hilbert filtering on the raw profiles around ${\it\omega}_{0}$ over $20T_{0}$ . The modulus of the velocity filtered for the numerical simulation with $a=1.6~\text{mm}$ is shown with a thick dashed line in figure 7(b). As the raw profiles, it is also very close to bidimensionality. By looking to the real part of the filtered signal, one can be convinced that the profiles for $v_{{\it\xi}}$ are exactly the same for positive and negative values. Hilbert transform allows us also to investigate the two-dimensionality of the TRI. As the numerical simulation with $a=1.6~\text{mm}$ is stable, one can use the one with $a=2.2~\text{mm}$ , which is unstable. Figure 7(b) shows the profiles along $y$ of the modules of $v_{{\it\xi}}$ filtered around ${\it\omega}_{0}$ (thick solid line), ${\it\omega}_{1}$ (thin solid line) and ${\it\omega}_{2}$ (thin dashed line) for this simulation. These profiles appear to be nearly constant in the $y$ direction, except close to the edges. This ensures that the instability which develops in the tank is also bidimensional.

To compare the numerics and the experiments, we can superimpose the experimental and numerical $v_{{\it\xi}}$ profiles after Hilbert filtering around ${\it\omega}_{0}$ for the different amplitudes of the wave maker. Figure 7(c) shows these different profiles: the solid blue line (respectively diamonds) corresponds to numerical (respectively experimental) data for $a=1.6~\text{mm}$ (respectively $a=1.5~\text{mm}$ ) while the dashed red line (respectively squares) corresponds to numerical (respectively experimental) data for $a=2.2~\text{mm}$ (respectively $a=3~\text{mm}$ ). Experimental points for $a=5~\text{mm}$ are represented by black circles. All data are normalized by the velocity in $y=0$  mm to avoid divergences in amplitude due to the location of the point where is measured the velocity along the branch and small discrepancies in wave maker amplitude. All data show again a very good two-dimensionality. The loss of amplitude on the edges, before the boundary layers is less than $10\,\%$ . Moreover, numerical and experimental data for $a=1.6$ and $a=1.5~\text{mm}$ (dashed blue line and diamonds) show a good agreement between each other. For higher amplitudes $a=2.2$ and $a=3~\text{mm}$ (red line and squares), where the TRI is developed, the agreement is also very good. Even for the experiment with $a=5~\text{mm}$ , where TRI is well developed, the shape of the experimentally measured transverse profile is in reasonable agreement with the simulation at $a=2.2$  mm. All numerical and experimental data collapsed in a reasonable way on a curve that ensures that all waves inside the tank are reasonably invariant in the transverse $y$ direction.

4.4 Generation of mean flows

We previously showed that the velocity field in the transverse direction for the wave attractor confirms the usual 2D approximation as a reasonable one. Nevertheless, it turns out that a closer look on the profiles of $v_{{\it\xi}}$ reveals a discrepancy between profiles of the maximum (max $v_{{\it\xi}}$ ) and the opposite of the minimum ( $-\!\min v_{{\it\xi}}$ ) of the velocity  $v_{{\it\xi}}$ as illustrated in figure 7(a). Indeed, the raw profiles for the attractor are not totally symmetric: when $v_{{\it\xi}}$ is positive (thick solid line), even if the profile looks 2D, there are some small peaks close to boundary layers which are not present in the profiles when $v_{{\it\xi}}$ is negative (thin solid line). This symmetry breaking suggests the presence of a mean flow, always in the same direction, close to the boundary layers and which modifies the velocity profiles close to boundary layers. To verify this, we analyse now further the data and check for possible mean flow generated within the tank. To extract this mean flow from the raw data, we used Hilbert filtering around ${\it\omega}=0~\text{rad}~\text{s}^{-1}$ .

The numerical simulation specially run to thoroughly study the mean flow has been performed with the same geometrical parameters as the one with the amplitude $a=1.6~\text{mm}$ for § 4.3. The attractor is stable and is precisely located at the same position in the tank as in figure 2. To limit the size of the data saved, the three components of the velocities were recorded only for the two horizontal planes, $z=100$  mm and $z=200$  mm, with 0.5 s timesteps. On these planes, the velocities were saved on a mesh with a good resolution of $128\times 256$ points.

Figure 8.

Mean flow (Hilbert filtering analysis at ${\it\omega}=0~\text{rad}~\text{s}^{-1}$ applied on numerical data) generated in the two different horizontal planes $z=200~\text{mm}$ (a) and $z=100$  mm (b). The colour represents the vertical velocity component, $v_{z}$ , while arrows represent the horizontal velocity components, $v_{x}\boldsymbol{x}+v_{y}\boldsymbol{y}$ , with the arrow scale displayed at the right of the figure. The wave maker is located at $x=0~\text{mm}$ , on the left of the figure, while the slope is indicated, on the right, by a vertical white line. The dashed rectangles delimit the regions where the attractor beams intersect with the planes located at $z=200~\text{mm}$ (a) and $z=100~\text{mm}$ (b).

Hilbert filtering at ${\it\omega}=0~\text{rad}~\text{s}^{-1}$ was applied on time series of the three velocity components for each point of the $x$ – $y$ plane. We consider the situation once the attractor reaches the steady state. Figure 8 shows the mean flow in the two different horizontal planes, $z=100$  mm and $z=200$  mm. The colour represents the vertical velocity component  $v_{z}$ , while arrows represent the horizontal velocity components, $v_{x}\boldsymbol{x}+v_{y}\boldsymbol{y}$ . The wave maker is located at $x=0$  mm, on the left of the figure. The slope is indicated, on the right, by a vertical white line. Figure 2 shows the location of the attractor in the $x{-}z$  plane. At $z=100$  mm, the attractor is reflecting on the slope while, at $z=200$  mm, it is reflecting on the wave maker. These are the two intense regions in figure 8, close to the slope at $z=100$  mm and close to the wave maker at $z=200$  mm. As reflections occur in these regions, it can be difficult to extract conclusions for the mean flow. Nevertheless, it is still possible to examine the mean flow in the regions where the attractor beams intersect with the planes located at $z=200$  mm and $z=100$  mm. These intersection regions are delimited by the dashed rectangles in figure 8. The flow in the branches of attractors shows a clear tridimensional behaviour, with recirculation zones where the horizontal components of the velocity vectors are forming jet-like currents close to the vertical walls of the tank to compensate for the mean flow in the central part of the wave beam. It is important to emphasize that the magnitude of the mean flow is much smaller than the velocity of the fluid due to the waves. Indeed, for this attractor, the typical amplitude of the velocity oscillations in the most energetic branch at the intersection with the plane $z=200$  mm is about $2.5~\text{mm}~\text{s}^{-1}$ while the typical velocities of the mean flow are about $0.1~\text{mm}~\text{s}^{-1}$ , thus less than $5\,\%$ of the velocity due to the waves.

4.5 Dissipation in the bulk and in the boundary layers

To further examine the three dimensionality of the flow, we determine the spatial distribution of the dissipation in the tank. We first compute the dissipation

where ${\it\nu}$ is the kinematic viscosity while the strain rate tensor is defined as $\unicode[STIX]{x1D626}_{i,j}=(\partial v_{i}/\partial x_{j}+\partial v_{j}/\partial x_{i})/2$ where $v_{i}$ and $x_{i}$ are respectively the components of the velocity field and of the position. Then, for a given $y$ -transversal coordinate, we average on a $20\times 20$  mm $^{2}$ square to define the locally averaged dissipation

which has also been averaged, after 42 periods after the start of the forcing, over three time periods. This quantity is shown as an inset in figure 9.

Figure 9.

Horizontally integrated dissipation $\tilde{{\it\varepsilon}}(y)$ (see the definition in (4.3)), normalized by the total dissipation $\tilde{{\it\varepsilon}}(W/2)$ . Panel (a) corresponds to the forcing amplitude $a=1.6~\text{mm}$ , while (b) corresponds to $a=2.4~\text{mm}$ . The solid blue line corresponds to a region located within the most energetic branch of the attractor, while the dashed black line corresponds to a region close to the centre of the tank outside the attractor. The inset present the locally averaged dissipation $\overline{{\it\varepsilon}}(y)$ (in $\text{mm}^{2}~\text{s}^{-3}$ ) for both forcing amplitudes within the most energetic branch.

However, as for example in Bonometti, Balachandar & Magnaudet (Reference Bonometti, Balachandar and Magnaudet2008), it is even more explicit to integrate this result horizontally along the transversal $y$ -coordinate as follows

in which $W$ is the total width of the tank.

We plot in figure 9 the horizontally integrated dissipation $\tilde{{\it\varepsilon}}(y)$ , normalized by the total dissipation $\tilde{{\it\varepsilon}}(W/2)$ within the parallelepiped 20 mm $\times$ 20 mm $\times$   $W$ , in two different locations: one in the most energetic branch after the focusing onto the slope (with the subscript $a$ for attractor), the second one around the centre of the tank (with the subscript $c$ for centre) and, thus, far from any branch of the attractor. For a weak forcing $a=1.6$  mm which leads to a stable attractor, one gets $\tilde{{\it\varepsilon}}_{a}(W/2)=1.7$ and $\tilde{{\it\varepsilon}}_{c}(W/2)=0.013~\text{mm}^{2}~\text{s}^{-3}$ . For a larger forcing $a=2.4$  mm leading to unstable regime, above values are typically multiplied by a factor two, since one gets $\tilde{{\it\varepsilon}}_{a}(W/2)=3.7$ and $\tilde{{\it\varepsilon}}_{c}(W/2)=0.028~\text{mm}^{2}~\text{s}^{-3}$ .

Above figure 9 emphasizes that, within the most energetic branch, the spatial distribution of the dissipation is only weakly changed while increasing the amplitude of the forcing, and therefore by passing from a stable to an unstable attractor. Approximately a quarter of the dissipation is located in the boundary layers, and the three remaining quarters are within the bulk of the flow. The dissipation outside the attractor is significantly altered by the amplitude of the forcing. This is presumably because TRI has generated secondary waves which propagate, outside the attractor and more precisely toward the centre of the tank. However, it is important not to overestimate this result since the total dissipation in the centre of the tank is more than 100 times smaller than its counterpart within the most energetic branch of the attractor. It was of course expected, since this is where waves, and therefore energy, are trapped. Moreover, figure 9 reveals that, in the most energetic branch of the attractor, approximately 25 % of the total dissipation occurs within the boundary layers which occupy less than 10 % of the total width.

4.6 Synthetic schlieren test

Using the numerically computed density gradients in the transversal direction, we can provide an estimation of the error made on experimentally measured density gradients by assuming that the flow is bidimensional inside the tank. We use the simulation for $a=2.2~\text{mm}$ as in § 4.3, by considering the density gradients instead of velocity fields. These data were available on a point located on the most energetic branch, for all $y$ . Figure 10(a) shows the time history of the horizontal density gradient in $y=0$  mm. TRI occurs after $35T_{0}$ in this simulation.

Dalziel et al. (Reference Dalziel, Hughes and Sutherland2000) have shown that the angle of deflection of a light ray crossing the tank in the $i$ direction is

$i$ being $x$ or $z$ . Here $L$ is the width of the tank and $n$ the optical index of water. The quantity $\text{d}{\it\rho}/\text{d}n$ is essentially constant and equal to $4.1\times 10^{3}~\text{kg}~\text{m}^{-3}$ . Equation (4.4) assumes that the flow is bidimensional and that the density gradients do not depend on the $y$ direction. This deflection induces displacement of the dots on the camera CCD from which we can measure the density gradients.

Figure 10.

(a) Time history of horizontal density gradient in $y=0~\text{mm}$ for the numerical simulation with $a=2.2~\text{mm}$ . (b) Displacement of the dots on the CCD camera that the simulation would have caused if the density gradients have been measured using the synthetic schlieren technique at $t=30T_{0}$ , before the onset of TRI. (c) Same as (b) after the onset of TRI ( $t=180T_{0}$ ). In both cases, the displacements are plotted for two periods of the wave attractor around the specified time. The blue circles indicated the displacement caused if the flow is assumed to be 2D, while the red squares show the displacement caused without this assumption.

With the numerical simulations, one can proceed in the inverse way. Knowing the density gradients and especially their dependence in $y$ , one can compute the displacement that these gradients would have caused on a CCD camera placed at typical experimental distances. For density gradients depending slightly on $y$ , one has to integrate over the width of the tank to get the deflection

This integration assumes that the deflection of the ray is small in $x$ and $z$ directions, so the ray remains more or less at $(x,z)$ constant. That is why one can use only a line profile in $y$ for density gradients. In typical experimental conditions, the CCD camera is placed at $1.75$  m of the tank and at $2.25$  m of the random dot pattern. The camera has a lens of $25$  mm. Considering these experimental parameters, simple geometrical calculations then allows us to convert the angle deflection ${\it\alpha}_{i}$ into the displacement on the camera ${\it\Delta}_{i}$ . Displacements are indicated on figure 10(b,c) by the red squares, before the onset of TRI (b) and after the onset of TRI (c). The order of magnitude of these displacements is of few micrometres, while the pixel size of the two different cameras used are equal to $3.45$ and $4.65~{\rm\mu}\text{m}$ . Thus, the displacements are of few pixels and sufficient for data treatment. Interestingly, one can also compute the displacements as if the flow was strictly bidimensional, by assuming that, once given $(x,z)$ , the density gradients in $y=0$ are the same for all $y$ in the tank. The deflection is thus as in (4.4). These displacements are indicated in figure 10(b,c) by the blue circles, before the onset of TRI (b) and after the onset of TRI (c). The difference between blue circles and red squares is very small in both situations, before and after the onset of TRI. This shows that the bidimensional assumption for the flow is robust and the error using it is very small. One can see in figure 10(b,c) that the displacements of blue circles are a slightly larger than those for red squares. This shows that the density gradients in the centre plane $xOz$ are slightly underestimated in the experimental measurements performed with the synthetic schlieren technique, typically by a factor of $5\,\%$ without TRI and by a factor of $7\,\%$ with TRI. This combined experimental–numerical check is a nice confirmation of the validity of the synthetic schlieren technique that justifies the comparisons performed in §§ 3.1 and 3.2 of the present paper.