Three dimensional flows beneath a thin layer of 2D turbulence induced by Faraday waves

Faraday waves occur on a fluid being subject to vertical shaking. Although it is well known that form and shape of the wave pattern depend on driving amplitude and frequency, only recent studies discovered the existence of a horizontal velocity field at the surface, called Faraday flow. This flow exhibits attributes of two-dimensional turbulence and is replicated in this study. Despite the increasing attention towards the inverse energy flux in the Faraday flow and other not strictly two-dimensional (2D) systems, little is known about the velocity fields developing beneath the fluid surface. In this study, planar velocity fields are measured by means of particle image velocimetry with high spatio-temporal resolution on the water surface and at different depths below it. A sudden drop in velocity and turbulent kinetic energy is observed at half a Faraday wavelength below the surface revealing that the surface flow is the main source of turbulent fluid motion. The flow structures below the surface comprise much larger spatial scales than those on the surface leading to very long-tailed temporal and spatial velocity (auto-) correlation functions. The three-dimensionality of the flow is estimated by the compressibility, which increases strongly with depth while the divergence changes its appearance from intermittent and single events to a large scale pattern resembling 2D cut-planes of convection rolls. Our findings demonstrate that the overall fluid flow beneath the surface is highly three-dimensional and that an inverse cascade and aspects of a confined 2D turbulence can coexist with a three-dimensional flow.


Introduction
Faraday waves (Faraday 1831) are subject to studies for a large variety of applications, ranging from bio-medicine to material sciences (e.g. controlled pattern formation, walking and orbiting of droplets) (Saylor and Kinard 2005;Couder et al. 2005). A complex and random transport of floating particles has been attributed to non-linear interactions at the surface of the Faraday wavefield, such as imperfections and traveling waves (Saylor and Kinard 2005;Ramshankar et al. 1990). However, only recent studies characterised the horizontal fluid flow produced by the waves, called Faraday flow. It was shown to exhibit attributes of two-dimensional (2D) turbulence, namely an inverse energy cascade (von Kameke et al. 2011(von Kameke et al. , 2013Francois et al. 2013). The presence of an inverse energy cascade in 2D turbulence has been theoretically predicted (Kraichnan 1971) and was observed in both numerical and experimental results for different fluid flows (von Kameke et al. 2011;Farazmand et al. 2011;Boffetta and Ecke 2012;Liao and Ouellette 2013;Francois et al. 2013) and references therein). Energy is introduced at intermediate forcing scales and transferred upwards to larger scales, resulting in a net inverse energy flux. Under particular conditions, this phenomenon can even lead to energy condensation, by which large and ordered flow structures emerge from the seemingly disordered motion at small scales (Xia et al. 2009;Musacchio and Boffetta 2019;Shats et al. 2014).
Due to the turbulent velocity field developing at the surface, the Faraday flow offers a test-bed to study the surface dispersion of objects with distinct geometries (Francois et al. 2018;Xia et al. 2019;Yang et al. 2019a, b), opening up a new study ground for fluid-structure and fluid-particle interactions. It was recently discovered in Francois et al. (2020) that the mean power extracted from a rotor placed in a strongly turbulent Faraday flow is highly dependent on the rotor thickness, influencing the energy of the angular velocity fluctuations. However, the details of the velocity field below the surface which could lead to this dependency were not investigated. Recent studies have been conducted to assess the three-dimensionality effects in quasi-2D flows, either electro-magnetically-driven (Kelley and Ouellette 2011;Martell et al. 2019), or produced by parametrically-excited waves Xia and Francois 2017). Complementarily, the occurrence of an inverse cascade in thick layers and the coexistence of 2D-and 3D-turbulence is also subject of current investigations (Biferale et al. 2017;Kokot et al. 2017;Biferale et al. 2012).
The details about the creation and the sub-surface velocity field of the Faraday flow, however, are to date not well understood. Especially intriguing is that the 2D turbulence at the liquid surface is produced by agitation in the vertical direction. It is likely that the Faraday waves induce this pronounced turbulent surface flow as indicated by data from Xia and Francois (2017) and recent theory that relates the creation of solenoidal currents on the fluid surface to its elevation for two plane surface waves (Filatov et al. 2016). For a thick layer of 3 cm liquid height the system is by no means shallow in comparison to the typical scale of energy injection (of the order of the half Faraday wavelength ≈ 5 mm) and, as we show in this study, pronounced three dimensional flows occur in the bulk. Here, we take an initial step in understanding the interplay of 2D turbulence at the fluid surface and the 3D fluid flow beneath for the Faraday flow. For the first time, we measure the sub-surface velocity fields by means of planar PIV with a high spatial and temporal resolution at multiple horizontal planes at different depths below the surface.

Water container and shaker set-up
Faraday waves are investigated in a circular container of acrylic glass (diameter 290 mm), similar to those used in pioneering studies (von Kameke et al. 2011(von Kameke et al. , 2013Francois et al. 2013) filled with distilled water at 21 ± 1 • C. A filling height of 30 mm is chosen for a deep water approximation, such that the depth is larger than twice the wavelength of the Faraday waves at the fluid surface. The independence of the results from the water filling height was assessed by performing additional test measurements with a filling height of 35 mm. These strongly resemble the results at 30 mm height, which confirms the deep-water approximation. The results at a filling height of 35 mm are reported in the supplementary materials, Figs. S5 and S6. The container is vertically shaken by an electromagnetic shaker (TIRA TV5220). A schematic representation of the experimental set-up is shown in Fig. 1.
Monochromatic sine forcing at f f = 50 Hz is imposed on the shaker from a function generator (RIGOL), which results in a Faraday wavelength F of 9.5 ± 1 mm and a subharmonic Faraday wave frequency response of 25 Hz. The acceleration of the container is measured with an accelerometer (Kistler, ±5 g, sensitivity 1 g/V±10% ) and read out by a high-frequency digitizer (Spectrum). Two levels of forcing acceleration a f are used: 0.47 and 0.70 g. The onset of Faraday waves is observed at a threshold acceleration of a th = 0.29 g, and determined as described in Ezersky et al. (1985) and Tufillaro et al. (1989). To account for hysteresis, the threshold acceleration was averaged using values observed by gradually increasing the forcing from standstill and those from reducing it from a stronger forcing. The resulting values for the dimensionless forcing amplitude, defined as 'supercriticality' = (a f − a th )∕a th , are = 0.61 and = 1.41 , respectively. In this study, measurements are carried out for a somewhat weaker forcing compared to similar experiments [where a supercriticality of = 1.7 was chosen , since this was found to be an optimal setting with regard to wave height and image quality for particle image velocimetry. Steep waves cause more reflections on the water surface and thus difficultate the particle imaging. The upper bound for the container acceleration in this experiments is the formation of water droplets forming on the surface of the waves (also called supercriticality) when = 3.48.

Camera and image acquisition
A second signal (standard TTL) from the function generator is used to trigger the high-speed camera (PCO dimax HS2, 12 bit depth). The camera is synchronised with the dominant frequency of the waves, which is found at the first subharmonic of the driving frequency f = f f ∕2 = 25 Hz, and records with a frame rate of 400 fps, corresponding to 1∕16 T = 0.0025 s, T being the Faraday wave period. The phase difference between the TTL driving the camera and the sine function fed into the shaker was then carefully monitored through the digitizer and tuned to capture the point of almost zero amplitude in the wave height (flat liquid surface) that occurs then at every eighth image. The camera is placed at the side of the shaker supports, and an optical prism-mirror is used to deflect the camera line of sight in the vertical direction. The camera resolution is 1400 × 1050 pixels, and images had to be cropped to 1400 × 850 pixels for the sub-surface measurements to match the width of the laser sheet. Images are saved in a 16 bit format (.b16), and subsequently converted back to a 12 bit format (bit depth of the camera).

PIV measurements
Two PIV techniques are used for the measurements at and below the surface respectively, which mainly differ in the choice of illumination light source and tracer particles employed. For the measurements beneath the waver level, red fluorescent polyethylene microspheres are used (diameter of 10-45 μ m, Cospheric), illuminated by a continuous wave argon-ion laser (wavelength of 488 nm, Ion Technologies). An optical arrangement is used to deflect the laser beam (first upwards and later again horizontally) to the desired measurement height. Afterwards, the light passes through light sheet optics (ILA 5150 GmbH) to generate a light sheet of 60 mm in width and 1 mm in thickness. The thickness of the laser sheet is much larger than the maximal vertical displacement of the container, which was analytically estimated and measured to be less than 0.0071 mm. Further, since the imaging is triggered to the container movement the measurement is always taken at the same height from the container bottom for one measurement height configuration. The tracer particles have a density of 0.995 g/cm 3 , and uniformly disperse in the water volume, when additionally treated with a surfactant, as described below. The particles have a Stokes number St ∼ O(10 −3 ) , which indicates that the inertial effects are negligible and particles can be considered as fluid tracers (Ouellette et al. 2008). The particles are neutrally buoyant since after a day of disposal no accumulation of tracers at the bottom or at the surface is observed. A high-precision long-pass filter (Edmund Optics) is used to capture the fluorescence of the particles (peak at 607 nm) and simultaneously shield the camera sensor from the laser light. The image of the fluorescent particles on the camera chip results to be 2-6 pixels in diameter. The described settings were used to measure the velocity fields at multiple horizontal planes with height h from the container floor, h = 3, 4, 5, 10, 15, 20, 25, 26, 28 mm for the weaker forcing at a f = 0.47 g, and h = 4, 10, 21, 27, 29 mm for the stronger forcing at a f = 0.70 g, see Fig. 1b). However, due to total light refraction at the water surface ( h = 30 mm), the combination of laser and fluorescent particles could not be used to measure the velocity field right at the water surface itself. In this case, a combination of floating hollow glass microspheres (diameter of approx. 70 μ m, St ∼ O(10 −3 ) , Fibre Glast) and a back-light (LED panel) was employed instead.  (4) is measured with accelerometers (5). A prism-mirror (6) is used to deflect the camera field of view. All the signals are monitored with a digitizer (7), data is saved on a laboratory computer. Numbers (8a) and (9a) depict the laser and its optics, (8b) shows the LED panel for the backlight PIV. b Reference for the height h of the measurement planes, measured from the container bottom For both PIV techniques, 0.3 g of particles are wetted in a 10%-solids solution with a surfactant (1% Tween 80 solution, Polysorbate 80, non-ionic surfactant). This helps to uniformly disperse the naturally buoyant particles (fluorescent) in the water volume, and the same surfactant is used in all measurements in order to avoid differences in the waves (e.g. due to changes in surface tension). With the available camera resolution (1400 × 1050 px), the conversion factor for the spatial calibration of the field of view (FOV) at the fluid surface h = 30 mm is 18.9436 px/mm for the case with a f = 0.47 g, resulting in a FOV of 73.90 × 55.42 mm 2 (details on the FOVs of the other settings are found in the supplementary materials, Tables S1 and S2). Figure 1 provides a schematic representation of the experimental set-up for the two PIV techniques described above, and Fig. 2 shows an example of raw images of the fluorescent particles at h = 29 mm, a f = 0.70 g. PIV measurements of the Faraday flow are carried out on the surface and at different depths in the water. Subsequently, the data is analysed using PIV View 2.6 with suitable interrogation window size and time intervals to have 4-6 particles per window and a particle displacement of 5-8 pixels. Due to the temporally well resolved measurements, these values could be adjusted by adapting the temporal step between successive images loaded to the PIV algorithm. The z-component of velocity (normal to the light sheet planes) could not be reconstructed from the available set-up. For the following figures and diagrams, the notation u = (u, v) ⊤ will be used to denote the velocity field and its components in x -and y -direction respectively, and h will be used for the height of the measurement plane with respect to the container bottom (see Fig. 1b). The background image in Fig. 3a) corresponds to the actual background-corrected images captured by means of backlight shadowgraphy, whereas in Fig. 3b the original background-corrected images have been inverted for better visualization of the velocity arrows. Measurement uncertainties and error estimation for the PIV techniques are quantified with regard to the following aspects. First, the uncertainty in manually setting the measurement height is estimated at h = ±0.5 mm. Further, calculation of particle displacement x PIV and the determination of the spatial conversion factor from the calibration target x cal are typical sources of error in the PIV algorithm.
The particle displacement error is estimated as explained in Raffel et al. (2014) and found to be x PIV ≤ 0.1 pixels. This in turn yields a maximum uncertainty in the u -and v -velocity estimation of u PIV = 0.16 mm/s for the surface measurements at higher forcing ( a f = 0.70 g). The error from the spatial calibration is estimated to be x cal ± 2 pixels over the total length of the calibration target, which ranges from 1114 to 1371 pixels. From this value a maximal uncertainty of u cal = 0.0011 mm/s can be estimated on the surface, which is however negligible compared to the error from the particle displacement calculation and therefore only u PIV is shown. The uncertainty estimation is included as error bars in Fig. 6a. Figure 3 shows an example of instantaneous velocity fields for the a f = 0.47 g case at the water surface ( h = 30 mm, Fig. 3a), and three heights below the surface: ( h = 28, 25, and 4 mm, Fig. 3b-d). The background images provide visual validation of the PIV calculations. The color code of the vectors indicates the absolute values of the velocity vectors � u � = √ u 2 + v 2 (also called magnitude of velocity, norm or vector length). Three distinct flow regimes can be identified from the evolution of the velocity-fields in Fig. 3a-d. The first regime is the Faraday flow on the water surface, with the characteristic presence of multiple vortices with variable length scales, as also observed in von Kameke et al. (2011), as well regions of jet-like flow in which the flow is strongly accelerated, similar to the riverlike structures defined as "trajectory bundles" in Francois et al. (2018). ere, an inverse cascade and a characteristic energy spectrum, the imprint of 2D turbulence can be detected as shown in Fig. S1. Second, a transition regime can be observed in a thin layer with a depth of F ∕2 below the surface where the structures gradually become slower and less turbulent but still prevail in its vortical appearance. These vortices seem to be tightly related to the surface and are not observed for h < 25 mm where the third regime is observed. Here, the velocity fields comprise larger flow structures oriented in dominant directions in the available field of view throughout the entire measurement time. The characteristic size of the flow pattern at different depths is quantified in Fig. 4, which depicts the time-averaged spatial correlation lengths L 0.5 for the u−component of the velocity. The 2D spatial correlation parameter u is calculated in Fourier space as presented in Fig. 2 Example of raw images of the fluorescent particles from the case h = 29 mm, a f = 0.70 g. A zoomed-in region of 350 × 130 pixels is shown. The particles appear with a diameter of 4-6 pixels on the camera chip. Faraday wavelength F depicted for reference Jähne (2005) and converted back to space domain, normalized with its maximum and time-averaged over all available time steps, here denoted by (⋅) . L 0.5 is finally determined by averaging the distances between the peak correlation at zero spatial shift, where u = 1 , and the points on the contour where u = 0.5 . It can be observed that in a thin layer with approximate thickness F underneath the water surface the correlation length is about half the Faraday wavelength, whereas at further depths the characteristic length of the structures consistently grows. The error bars in Fig. 4 show the standard deviation of L 0.5 along the contours for u = 0.5 , and large deviations can be observed for h < 25 mm. This strong correlation for a distinct spatial direction is due to the previously mentioned dominant direction of the flow below the surface throughout each experiment. For the higher amplitude forcing ( a f = 0.70 g) at h = 10 mm the correlation lengths were larger than the size of the field of view for all measurements. Therefore they could not be determined within this measurement campaign but would A dramatic velocity reduction can be appreciated below the surface. In fact, the mean absolute velocity drops to 1∕e of the surface value in 1.5 mm and 3 mm for a f = 0.47 g and a f = 0.70 g, respectively. At half the Faraday wavelength (5 mm underneath the surface) the mean absolute velocity has dropped by a factor of about 6.5. For values of h below 20 mm the mean absolute velocities level out at a nearly constant plateau. The trend in the profiles shown in Fig. 5a could be well fitted with an exponential function, ⟨� (h)�⟩ = a + b ⋅ exp(c ⋅ h) , using the built-in Matlab curve fitting tool. The quality of the fit was analysed by the coefficient of determination R 2 ≥ 0.98 , which is the square of the correlation between the measured values of ⟨� �⟩ and the predicted values from the fit. A value of R 2 = 1 implies a perfect fit, while a value R 2 = 0 means that the fit can not explain the variation of the measured values with height (Devore 2011). The decay rates of the exponential fits are c = 0.36 mm −1 for a f = 0.47 g and c = 0.76 mm −1 for a f = 0.70 g.

Velocity fields
The exponential decay of velocity magnitude is a somewhat unexpected result, since it corresponds to the most classical results of horizontal velocities beneath travelling water waves (Dietrich et al. 1980;Breivik et al. 2014) and the authors are not aware of a theoretical explanation for this observance in the Faraday flow. Moreover, the decay rates of the velocity norm are of the same order of magnitude of the exponential decay found in the theory of traveling waves, namely the wavenumber k F = 2 ∕ F ≈ 0.66 mm −1 . In Fig. 5a, the error bars indicate the standard deviation from the mean of the absolute velocities, computed over all time steps and grid points, which gives an idea about how much larger the turbulent fluctuations are on the water surface compared to the sub-surface flow fields (note the semilog scale) . This is confirmed by analysing the turbulent kinetic energy ( TKE ), computed as TKE = 1 2 u �2 RMS + v �2 RMS , where u ′ RMS and v ′ RMS are the RMS values of the velocity fluctuations u � = u − ⟨u⟩ , presented in inset graph Fig. 5b. The flow structures developing beneath the turbulent surface layer barely show any turbulent kinetic energy below a depth of approximately half the Faraday wavelength. The decrease in TKE is not merely caused by the decrease in mean flow velocity, but by the presence of large-scale and slower structures that persist longer in time as can be visually observed in the supplementary movies M1-M3 and is further analysed in Fig. 7.   Figure 6 depicts the profiles of the root-mean-square (RMS) values of the velocity components, u RMS and v RMS , for both forcing amplitudes. These follow the trend of exponential decay as evaluated for the mean of the absolute velocities. Furthermore, they reveal that the differences of RMS velocity values for u RMS and v RMS increases below the surface. This was further investigated by looking at the probability distributions of the velocity components. The inset graphs in Fig. 6b, c show the combined probability distribution of u and v for both forcing amplitudes and at two different measurement heights ( h = 30 mm and h = 10 mm respectively). The distributions are averaged through all available timesteps and grid points ( ≈ 10,000 ) for the surface measurements, here u -and v-components are combined together to improve the statistics). The symmetric nature of the velocity fields on the surface can be appreciated, since for both forcing amplitudes the velocity probability shows a symmetric distribution with respect to the origin (Fig. 6b). The curve for the stronger forcing is also somewhat flatter since a broader range of velocities is reached and the TKE is higher. The situation is considerably different at h = 10 mm, where for the weaker forcing (blue bars) two distinct peaks can be recognized, which indicate different mean components of the velocity in u and v , and a broad curve is seen for the stronger forcing with non-zero mean. For both forcing amplitudes the velocity probability distributions are asymmetric over the available field of view and measurement time-span in horizontal planes below the surface. This indicates that the temporal persistence of the structures surpasses the experiment time, which is further analysed by computing the velocity autocorrelation function, see Fig. 7 and the next section. It needs to be mentioned that the average velocity in a bounded flow is expected to be zero if computed over the entire domain size or a much higher number of statistically independent realizations of the same experiment, but both surpasses the possibilities of the current experimental set-up and data storage capabilities.

Velocity autocorrelation function
To quantitatively analyse the long temporal persistence of the flow structures observed by the eye (in Fig.3c, d and also supplementary movies M1-M3), the autocorrelation function of the velocity signals at different heights is analysed. The velocity autocorrelation function is calculated in Fourier space as presented in Jähne (2005) and converted back to the space domain. Spatial averaging (denoted with ⟨⋅⟩ ) is carried out over all the grid points and every curve is then normalized with its maximum (at 0 time lag) to obtain the averaged autocorrelation coefficient ⟨ ⟩ . The results are presented in Fig. 7 for the u-velocity fields, namely ⟨ uu ⟩ , at different values of h for both forcing amplitudes.
Very similar results have been obtained for the v-velocities (not shown here). In most cases at a f = 0.70 g, depicted in red color scale and dashed line, data could be gathered over 2-3 replications of the experiments, which were statistically fully independent (the experiment has been restarted). To further improve the statistics for the surface case the autocorrelation function of u and v are averaged, as the velocity field is symmetric with respect to the origin, as shown in Fig. 6b. From the autocorrelation function diagram a trend can be immediately recognized, namely that the velocity structures change on short time scales in the proximity of the surface, and that the decorrelation of velocity signals occurs much faster for stronger forcing, dropping below the 1∕e-value in under 2 s for the a f = 0.47 g case, and in under 1 s for a f = 0.70 g. The velocity fields, however, become extremely time-persistent close to the container bottom (lower values of h ) for both forcing amplitudes. This trend is found to be consistent with the results from the velocity profiles from Figs. 5 and 6: at the water surface the flow is faster, with stronger turbulent fluctuations and structures that decorrelate more quickly in time. However, below a certain height h (25 mm for a f = 0.47 g and 21 mm for a f = 0.70 g) the high values of correlation indicate the presence of slower structures with time scales much longer than the time of the recordings.
In fact, for the lower forcing at a f = 0.47 g the distribution of the wavefield close to the borders shows stripy patterns oriented in the radial direction, which in turn could lead to the development of different velocity layers below the surface. However, by increasing the forcing to a f = 0.70 g, a highly disordered isotropic wavefield is achieved on the whole surface, while the strong temporal velocity autocorrelation at low heights persist through the analysed timespan. This confirms that the large-scale structures below the surface develop independently from boundary effects in the wave pattern. It can also be noted that the autocorrelation function at the surface for the stronger forcing presents an oscillating behaviour, which might be related to the stronger vertical motion of the waves. This behaviour is difficult to interpret by looking at velocity information on horizontal planes. It is suspected that a small portion of the horizontal motion of the particles on the surface, is due to the wave motion, since at high forcing amplitudes the surface is never entirely flat even with a very careful synchronization of the camera recording to the vertical agitation. Figure 7b shows the decorrelation time of the velocity autocorrelation function over height h extracted from Fig. 7a by taking the time at which ⟨ uu ⟩ = 0.5 . For both forcing amplitudes the decorrelation time gradually increases below the surface. The profile is not easily described by an inverse proportionality to 1∕u which would be expected for a mere decrease due to slower flow structures. Rather, the steep increase for lower measurement heights indicates that the flow structures additionally become more persistent in time. In fact, for h < 20 mm the correlation time could not be quantified with the current setup since the camera buffer memory only allowed a maximal measurement time of 12 s at the necessary frame rate.

Vorticity and divergence fields
The relative size and behaviour of the spatial structures of the velocity fields is further investigated by analysing instantaneous contours of vorticity and divergence, computed from the two-dimensional velocity data as z = v∕ x − u∕ y and divu = u∕ x + v∕ y respectively. The results are presented in Fig. 8 for the a f = 0.47 g case and in Fig. 9 for a f = 0.70 g. For the higher forcing, the velocity fields have been slightly smoothed with a Gaussian filter (with a kernel of 3 × 3 velocity grid nodes) to remove measurement noise prior to calculating the sensitive gradients. The results reflect the findings from the velocity profiles and velocity autocorrelation functions. On the water surface, shown in Figs. 8a and 9a, regions of alternating vorticity peaks are densely distributed across the entire field of view, with a size varying mostly in between F ∕2 and 2 F . Below the surface, the vorticity intensity drops significantly by more than one order of magnitude for both forcing amplitudes, as shown in Figs. 8c and 9c. The vorticity structures become larger and smoother and no obvious imprint of the Faraday waves is observed anymore.
The divergence fields depicted in Figs. 8 and 9b, d draw a different picture. The divergence on the surface is overall lower than the vorticity, especially for the lower forcing case, further it presents only localized peaks and spatial scales much smaller than that of the vorticity structures. The divergence fields for subsurface planes at h = 4 mm in Figs. 8 and 9d indicate the presence of persistent vertical motion that, in the horizontal-plane projection, appears as a lattice of sinks and sources and resembles the pattern of a 2D cut-plane through convection rolls. The sinks are in this case narrow stripes of strong negative divergence, whereas the regions of positive divergence (flow from the bottom to the surface) are weaker but larger in size, and present an elliptical shape with a minor axis of approximately 2 F . This divergence pattern underlines the strong three-dimensionality of the flow that, to our knowledge, has not been previously reported. The only quantitative work on flow patterns and velocity fields below the surface of the Faraday waves that we are aware of considers only purely longitudinal Faraday waves confined in a narrow channel and for shallow layers (Gutierrez et al. 2016). The reported patterns resemble the ones observed in this study qualitatively with the important difference that in our study the stripes additionally exhibit a temporal evolution and are disordered in space. Despite the non-zero value of divergence at the surface of the Faraday flow, it was shown in Xia and Francois (2017) that the divergence (quantified by the compressibility parameter) rapidly decorrelates over four Faraday wave periods. This implies that the divergence does not, in the mean, contribute to the particle displacement of surface bound particles and, in turn, confirms the quasi-two-dimensional nature of the Faraday surface flow even for thick layers. The compressibility parameter has been analysed at different depths, as shown in Fig. 10. As presented in Xia and Francois (2017), the compressibility parameter is computed as: where ⟨⋅⟩ T av denotes the averaging over the time, where T av = n ⋅ T is the number of Faraday wave periods. Figure 10  Fig. S3 in the supplementary materials. For both forcing amplitudes the compressibility increases below the water surface, which is in accordance with the results seen in Figs. 8 and 9, where the relative importance of the divergence over the vorticity increases with depth. The important mark of a compressibility ⟨C⟩ = 0.5 that corresponds to the value for a 2D flow entirely made up of sinks and sources ) is crossed at a depth of approximately one Faraday wavelength. This shows that the flow becomes very three-dimensional for larger depth.

Conclusions and outlook
The measurement and analysis of the Faraday flow at and below the water surface have revealed a rich variety of complex flow structures and unveiled the existence of previously unknown three dimensional bulk flows. This finding immediately raises the question about the mechanism fueling the fluid motion: does the bulk flow drive the surface flow or does the surface flow drive the bulk flow? To address this question, the turbulent kinetic energy, the mean of the absolute values of the velocities and the velocity RMS were studied at different horizontal planes beneath the water surface. An impressive exponential decrease with depth was observed for these measures and also for the vorticity RMS (see supplements), proving that the forcing mechanism is confined to a thin layer of approximately half a Faraday wavelength below the surface. This confinement to the fluid surface has currently been suggested by theoretical and experimental work on wave-driven flow Filatov et al. 2016;Levchenko et al. 2017;Francois et al. 2015Francois et al. , 2017. Notwithstanding, the full picture of the 3D flow structures is complex. The flow structures on the surface are symmetric and their timescales and lengthscales are short, whereas below a surface layer of half a Faraday wavelength large flow structures are observed that only vary on long timescales. From our analysis, we can differentiate three different vertical regimes, namely the 2D Faraday flow at the surface, the transitory regime of decay of vorticity and the regime of large and persistent structures with dominant divergence patterns. Our findings in this regard shed light on the emergence of 2D turbulence and its imprints such as an inverse cascade in three-dimensional flows. A subject of that has attracted increased attention recently (Biferale et al. 2012) with applications in the fairly new fields like the generation of active turbulence (Kokot et al. 2017). It is currently unknown which mechanisms causes dissipation and prevent energy condensation, a pile-up of energy in the largest scale, in the current set-up. Due to the thick layer and the confinement of the high fluid velocities and the turbulent kinetic energy to the fluid surface, the classical bottom drag effect should be negligible. Furthermore, the increase in three-dimensionality of the flow as indicated by an increase of the compressibility with depth suggests that there might be a transport of energy to lower depths. This energy might then be dissipated by 3D mechanisms like a direct energy cascade. Observances by eye, when neutrally buoyant fluorescent particles are used, suggest that the high peaks of the divergence at the fluid surface stem from small and fast vertical jets that are pushed from the surface downwards and might actually drive the bulk flow. A future study will consider the full three-dimensional velocity field to further unravel the interaction of the Faraday waves with the 2D turbulence at the surface and the persistent 3D bulk flow.