Gas flow dynamics over a plunging breaking wave prior to impact on a vertical wall

We present an experimental study on the gas flow field development over a plunging breaking wave prior to impact on a vertical wall. The variability of wave impact pressure over repeated measurements is well known (Bagnold, 1939). The formation of instabilities on the wave crest are postulated to be the main source of impact pressure variability (Dias and Ghidaglia, 2018). However, the mechanism that results in wave impact pressure variability and the influence of the gas phase in particular are relatively unknown. The velocity field of the gas phase is measured with particle image velocimetry, while simultaneously the local free surface is determined with a stereo-planar laser-induced fluorescence technique. The bulk velocity between the wave crest tip and the impact wall deviates from the mass conservation estimate based on the velocity profile between the wave crest and the impact wall. This is caused by a significant increase of the local gas velocity near the wave crest tip. The non-uniformities in the seeding concentration accumulate near the wave crest tip and reduce the accuracy of the velocity measurements. However, the bulk velocity estimate is significantly improved with a fit of the velocity profile that is based on a potential flow over a bluff body. Additionally, the development of vortex is observed and quantified for two typical measurements with either a disturbance on the wave crest or a smooth wave crest. The circulation development is comparable to the formation and separation of a vortex ring, which results in a saturated vortex that separates from the wave crest (Gharib et al., 1998). Furthermore, the impact location of the wave tip is altered by the formation of secondary vortices. The secondary vortex enhances the lift locally and alters the trajectory of the wave crest tip, which may result in additional variability of the wave impact pressure. © 2021 The Author(s). Published by ElsevierMasson SAS. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).


Introduction
Stricter emission regulations result in new use cases of liquefied natural gas (LNG), which can serve as a ''transition fuel'' until greener options become economically viable. This growth of the LNG market brings about an increased demand for floating liquefaction facilities, storage facilities, and shipping solutions [1]. Additional challenges arise with the widespread use of LNG, such as the requirement for lower filling levels and increased containment capacity [2]. The filling levels in LNG containment systems are currently restricted to limit the occurrence of sloshing. Sloshing increases the number of extreme wave impact events that can potentially result in structural damage [3,4]. The development of containment systems is limited by the scaling involved in review on slamming [22]. The wave impact process is subdivided into elementary loading processes (ELPs), such as the direct impact, the jet deflection, and the compression of the entrapped or escaping gas [23]. Each combination of elementary loading processes results in a different wave impact pressure distribution, and wave impact shape (i.e., type). The wave impact types are either defined as a slosh, a flip-through, a gas pocket, or an aerated type of wave impact [4,8]. For example, the flip-through wave impact has been studied in detail. A flip-through impact shows significant vertical accelerations of up to 1500g [24]. A gas pocket type impact is often created by a plunging breaking wave. The horizontal impact of a plunging breaking wave often encloses a gas pocket. After wave impact the gas pocket decreases in size and starts to oscillate [25,26]. However, the scaling of these physical phenomena, such as wave impact pressures and wave impact type, from small-scale to large-scale is not trivial [18,21].
Dynamic similarity is generally not easily obtained during wave impact events when scaling from small-scale to largescale [21]. The relevant similarity parameters change when the wave approaches the impact wall, which warrants the use of ELPs to identify the similarity parameters at each step of the wave impact process [23]. For example, the global flow, where the wave is not influenced by the presence of the impact wall (i.e., the increase in pressure in the gas pocket and increase in flow from the enclosed gas pocket), can be Froude scaled [22,27]. On the other hand, the scaling of the local flow is not straightforward. The local flow is influenced by the strong gas flow over the wave crest and by the surface tension of the gas-liquid interface [28], the gas-liquid density ratio [29,30], the compressibility of the gas (i.e., the speed-of-sound) [27], the possibility of phase change [6,31,32], and the aeration of the liquid [10,33,34]. The scaling of the local flow is not well understood, but especially the formation of ligaments and droplets are thought to be relevant for the variability in wave impact pressures [22,35].
The source of impact pressure variability in repeated wave impact experiments is thought to be the instability development on the wave crest. The mechanisms that are responsible for the formation of these instabilities are still largely unknown [35]. The shear force induced by the gas flow can result in the variability of the local flow [30,[35][36][37]. Additionally, the wave tip of the plunging breaking wave can be deflected [38,39]. Prior to impact, gas cushioning (i.e., the increase in pressure in front of the wave tip) can also result in deformation of the wave tip [40][41][42]. The wave tip deflection is shown to depend on the density ratio and the scale of the experiment [29,30]. In a previous manuscript the variability of repeated wave impact experiments was accurately determined, and a mechanism for free surface instability development was presented [43]. However, the influence of the gas flow on the local flow during repeated wave impact experiments has not investigated yet.
Developing waves generate a flow of gas over their free surface by moving through the stagnant air (e.g., similar to an obstacle moving through air). The wave crest geometry (i.e., the height and slope) and wave celerity define the gas velocity and possible flow separation [44]. In most measurements, flow separation occurs at the wave crest (e.g., the top of the wave), which results in a vortex [45]. The vortex often separates from the wave crest and the separated vortex lingers in the air [44]. Plunging breaking waves also develop a vortex over their wave crest as they plunge towards a free surface [46]. Additionally, the gas velocity at the wave crest tip is enhanced by the flow from the enclosed gas pocket. However, such experiments are often performed on small amplitude waves.
In the present study, particle image velocimetry (PIV) measurements are combined with accurate free surface measurements at both a local and global scale to further investigate the mechanism of impact pressure variability in repeated wave impact experiments. The wave generation method is repeatable on a global level, whereas the local flow displays significant variation in both the details of the wave crest shape and the impact location [43]. The local flow is investigated with both flow field (PIV) and stereo-planar laser-induced fluorescence (stereo-PLIF) free surface profile measurements. The PIV measurements are performed with two field-of-views (FOVs) to obtain both the global and local flow fields with sufficient accuracy. The velocity profile between the wave crest tip and the impact wall is repeatable. However, in Section 3.2 we will show that there is a mismatch between the bulk velocity based on the velocity profile and the mass conservation estimate. A velocity profile fit that resembles the potential flow over a bluff body is required to accurately determine the bulk velocity. Additionally, two typical vorticity fields are identified that depend on the shape and smoothness of the wave crest. Interestingly, the difference in wave crest shape also result in a change of the gas velocity field prior to impact. The change of the gas velocity field (i.e., the formation of a secondary vortex) results in additional lift at the wave crest tip. The formation of a secondary vortex near the wave crest tip is identified as a possible additional source of wave impact pressure variability.
This paper is organized as follows: the experimental setup and equipment are introduced in Section 2. The PIV method is described in detail, whereas references to the other elements of the experimental set-up are presented. Thereafter, the results are introduced and discussed (Section 3). First, in Section 3.1, the global gas flow around the wave crest is discussed for a typical measurement. Then, in Section 3.2, the velocity profile between the wave crest tip and the impact wall is discussed. Thereafter, in Section 3.3, the development of a vortex over the wave crest and near the wave crest tip are discussed in the context of wave impact location and possibly wave impact pressure variability. Finally, the findings are summarized in Section 4.

Wave flume
The measurements are performed in the wave flume of the Hydraulic Engineering Laboratory at the Delft University of Technology (Fig. 1). The flume has a width of 0.79 m and a height of 1 m with a cross-section of 0.79 m 2 . The total length of the flume is 39 m and the water depth is maintained at h 0 = 500.0 ± 0.5 mm for all measurements. The flume is equipped with a piston-type wavemaker that has a maximum stroke of 2 m. The time between measurements is reduced with an active reflection compensation system, that is only enabled after the impact of the focused wave on the vertical impact wall.
In this work, a plunging breaker with a large gas pocket that impacts on a vertical wall is studied. The large gas pocket wave is generated with a wave focusing technique [13]. The wave group components with their own phase velocity and wave length are focused on the focal point (x f ). The various wave periods of the waves in the wave group are clearly visible in the measured free surface elevation profile, that is shown in the inset of Fig. 1. The angle of the wavefront with respect to the impact wall and the impact type are determined by the focal point of the wave groups. A shift of the focal point results in either: a slosh, a flipthrough, a gas pocket, or an aerated impact [8]. A wave shape is experimentally determined that results in a high impact pressure while enclosing a large gas pocket. The nearly parallel front of the wave crest (i.e. a wave crest aligned with the impact wall) results in a high impact pressure [34]. The required wave shape is obtained with a trial-and-error approach, and a normalized focal The PIV camera measures the global wave shape and the velocity around the wave crest of the plunging breaker at the center plane of the flume. The velocity field around the wave crest is determined at two levels of magnification, respectively FOV1 and FOV2. The stereo-PLIF system measures wave crest details in a smaller field-of-view at the light-sheet location.
The current method of wave generation is similar to the large scale tests of the Sloshel project [18]. A similar depth-based Froude number is required to compare the measurements. The effective flume length is therefore scaled as λ = h 0 /h * 0 with h 0 the current water depth and h * 0 the reference water depth [23]. The reduced effective flume length is obtained by mounting a 20 mm thick transparent perspex wall on a rigid construction at a distance of L = 23.4 m from the wave board ( Fig. 1). However, exact Froude scaling is not achieved due to practical limitations (e.g., the camera measurement system limits the water depth to 500 mm). The Froude scaled ratio is (1 : 7.3) compared to the (1 : 6) ratio of the Sloshel experiments, which will result in a smaller wave (i.e., with a smaller gas pocket and lower wave impact height) [23,27].
The generation of nominal identical waves with a focusing technique is not trivial, as changes in the initial conditions, such as the water depth result in a different impact type [7,14,47]. The variance in impact type results from two sources of variability: the system variability, and the hydrodynamic variability. Minimization of the system variability is essential to study the hydrodynamic variability (e.g., the growth of free surface instabilities on a wave crest). The control parameter variation of the facility (i.e., system variability), such as the water depth, the piston motion, and residual motion are minimized within the limitations of the experimental facility. The wave generation is shown to be repeatable in a previous manuscript [43]. However, the comparison between measurements over several days is limited due to inevitable day-to-day variations present in wave flume facilities [14,27].

Experimental methods
The gas velocity and local free-surface profile are measured simultaneously with planar particle image velocimetry (PIV) and stereo planar laser induced fluorescence (stereo-PLIF). The stereo-PLIF technique is described in detail by van Meerkerk et al. [48], and the application of the technique in the current facility is described in [43]. In this manuscript the equipment used for the stereo-PLIF measurements is only introduced, and the interested reader is referred to the previous manuscripts for details concerning the data processing of the stereo-PLIF images [43,48]. The PIV equipment and processing are introduced and discussed in the following paragraphs.
Illumination of the PIV domain is provided by a high-speed Nd:YLF laser (LDY 304 PIV laser, Litron) with an output of approximately 10 mJ per pulse at a repetition rate of 2.5 kHz. A vertical light sheet is created at the center plane of the flume (Fig. 1). The light sheet is created by expanding a laser beam with a cylindrical lens of focal length −12.5 mm. Additionally, the beam is focused into a narrow sheet with a cylindrical lens (f = 1000 mm), which results in a light sheet thickness of approximately 2.5 mm.
The laser illuminates both the water droplets in the gas (see below) and the incoming wave. The water in the flume contains a low concentration of fluorescent dye (Rhodamine WT, Sigma-Aldrich at 120 mg m −3 ). The stereo-PLIF cameras are equipped with an optical high-pass filter that only shows the illuminated free surface. The free surface detected with the stereo-PLIF technique is projected onto the PIV plane. The PIV data is masked with the projected stereo-PLIF data in the final step of the PIV processing.
Seeding of the stagnant air in this experiment is not trivial. Small droplets or particles would precipitate onto the free water surface where they would affect the surface tension properties [44,49,50]. The effect of seeding on the gas-liquid surface tension is reduced with water droplet seeding [50]. A continuous supply of seeding particles is generated with a commercial ultrasonic humidifier (UHW, Medisana). The device is enabled prior to the experiments to fill up the flume with a fog of small water droplets. The diameter of the droplets (d p ) is determined with a interferometric particle imaging (IPI) technique and is approximately d p ≈ 10 µm (see Appendix A for details). The particle response time is approximately τ p = ρ p d 2 p /18ν f ρ f = 0.3 ms based on a particle density of ρ p = 998 kg m −3 , the density of air ρ f = 1.23 kg m −3 at standard atmospheric conditions, and the kinematic viscosity of air of ν f = 15 × 10 −6 m 2 s −1 . The settling velocity of the particles (e.g., ∆u = (ρ − 1)/ρgτ p ≈= 3 mm s −1 ) is small with respect to typical velocities of the wave impact. The tracer particles are estimated to follow velocity fluctuations up to a cut-off frequency of approximately 550 Hz [51,52]. Furthermore, the particles are able to adjust to the acceleration expected in the measurements based on the bulk velocity estimate [43,50]. The images are acquired with three high-speed CMOS cameras (Imager HS 4M, LaVision). The stereo-PLIF cameras (i.e., camera 1 and 2 in Fig for FOV2 (Fig. 2). The recording of FOV1 and FOV2 did not occur simultaneously, but rather serves as a highlight of differences in local gas flow velocity between globally repeatable waves. The frame rate (f aq ) is 2.5 kHz with an exposure time (∆t e ) of 363 µs.
The frame rate is obtained at a reduced image resolution off 1157 × 631 pixels for FOV1 and 1042 × 849 pixels for FOV2.
The PIV camera is calibrated for FOV1 with a PMMA plate with crosses equally spaced at a distance of 50 mm, whereas for FOV2 a two-plane dot-pattern target (Type 22, LaVision) is used. A third order polynomial mapping function is determined for both field of views to back-project the images to the center plane of the flume [52,53]. The error on the mapping function is 0.1 and 0.3 pixels at a scale of 2.1 and 6.3 pixels per millimeter for respectively FOV1 and FOV2.
The images are processed with a PIV algorithm implemented in Matlab (2019B, Matlab). The images are preprocessed with a min-max filter that normalizes the image contrast [52]. However, significant variation in image intensity remains due to both the reflection of light from the free-surface and the nonhomogeneous seeding of the gas. The variation in image intensity contaminates the correlation peak [54,55]. Therefore, symmetric phase only filtering [54] is applied with a symmetric filter to sharpen the correlation peak (see Appendix B). The images are processed with a multi-pass and multi-grid interrogation strategy with window deformation. The advantage of the multi-grid approach is the increased dynamic range of the PIV measurements [52]. The final multi-pass step results in an interrogation window of 16 × 16 and 8 × 8 at an overlap of 50% for respectively FOV1 and FOV2. The vector spacing is approximately 3.8 and 0.6 mm for the respective field of view's. Outliers are identified with a median test (typically 3% of the vectors in the initial pass) and replaced by interpolation [56].

Results and discussion
In this section the results are introduced and discussed. First, the global gas flow behavior over a plunging breaking wave prior to impact on a vertical wall is discussed and relevant definitions are introduced. Second, the gas velocity between the wave crest and the impact wall is determined, and the velocity profile is estimated from the valid velocity vectors. Third, the vorticity and circulation over the wave crest is determined. The vorticity field over the wave crest is discussed by highlighting two typical waves that show distinct behavior on both a global and local scale of the velocity field in the gas phase. Finally, the variability of two typical wave impacts is studied and a mechanism that enhances the variation is discussed.

Global gas flow
The development of the gas flow over a plunging breaking wave prior to impact on a vertical wall is shown in Fig. 3. This figure shows measurement M225, also defined as the perturbed wave, one of the two highlighted measurements in this work.
These measurements are part of a larger set of N = 21 measurements that were collected over two successive days. The panels show the stereo-PLIF data (Fig. 3d) and gas velocity fields superimposed on the original image of the wave with tracer particles. The gas is expelled from the gas pocket with crosssectional area A(t) and flows through the throat (∆x(t)) between the wave crest tip (⃗ x wt (t)) and the impact wall (Fig. 3a).
For t = −24.0 ms, the velocity from the gas pocket is mostly vertically directed (Fig. 3a). The flow structure resembles that of a wall-bounded jet, which is often observed in numerical simulations of wave impacts with a large gas pocket [57,58]. On the other hand, the mist of water droplets already shows some form of circulation near the tip of the disturbance (Fig. 3a). The disturbance on the wave crest acts as a roughness element that triggers flow separation and the development of a shear-layer. Unfortunately, the resolution of the FOV1 PIV data is not sufficient with a vector spacing of 3.8 mm to directly resolve the flow separation near the wave crest disturbance.
At t = −18.0 ms, the resolution of the PIV data is sufficient to detect a vortex at the back-side of the wave crest (Fig. 3b). A vortex is consistently formed in all wave impact measurements. The size of the vortex increases over time from panel b to d. Additionally, the gas velocity in the wall-bounded jet increases shortly before impact of the wave on the vertical wall. The development of the vortex resembles that of a plunging breaking wave impacting on a free surface [46]. However, here the plunging breaking wave impact on the vertical wall lacks the development of large-scale secondary vortices that were observed by Techet and McDonald [46] during the impact of a plunging breaker on a liquid free surface.

Local gas flow
The local gas flow is defined as the flow through the throat (∆x), that is the gap between the wave tip (⃗ x wt ) and the impact wall. The plunging breaking wave forces gas through the throat by reducing the cross-sectional area of the gas pocket (A(t)). The rate of change of the cross-sectional area (A) is assumed to be constant over the span-wise direction of the wave, whereas the throat distance (∆x) is expected to vary. Typically, the global flow of the wave is repeatable when the system variability is minimal. The repeatability of the global flow in this facility has been confirmed previously [43]. The power-law fit of the cross-sectional area was found to be A(t)−A(t = 0) = 1.74|t| 1.52 , which results in a rate of change of dA/dt = 2.64|t| 0.52 [43]. The throat size is determined by the velocity of the wave crest tip, which is here ∆x = 2.67t +C with C being a constant. The averaged or bulk velocity of the gas in the throat is, with the assumption of constant gas density, defined by mass conservation, as with dA/dt the rate of change of the cross-sectional area, ∆x the throat size, and t the time to impact [43]. A time shift (τ 0 ) is included in the bulk velocity derivation of the gas flow, as the exact moment of impact is difficult to determine from the stereo-PLIF and visualization images of the wave impact.
The bulk velocity of the gas flow in the throat is determined for all PIV measurements of both FOV1 and FOV2. The bulk velocity is obtained by integrating the velocity profile over three The field of view (FOV) of each measurement system is defined. The PIV measurements are performed at two magnifications, which results in a large field of view (FOV1) and a small field of view (FOV2) measurement. The stereo-PLIF (SPLIF) measurements are combined with the ellipse fit of the gas-pocket cross-sectional area [43] to mask the PIV images. The stereo-PLIF images of camera 1 and 2 (see Fig. 1) are combined into a single profile that describes the interface location over a domain of 175 mm × 120 mm. interrogation windows centered at the y-coordinate of the wave tip. The bulk velocity of the gas flow is averaged per FOV over all repeated measurements, which requires the x-coordinates of the wave tip to be matched. The difference between the x-coordinate of the wave tip over repeated measurements is minimal, and the maximum temporal shift over all repeated measurements is only 2.5 ms. The bulk velocity of the gas flow over time is shown for both FOV1 and FOV2 in Fig. 4. The bulk velocity is repeatable for both FOV1 and FOV2 with an averaged standard deviation over bulk velocity of 2.2% for both measurements.
The bulk velocity of the gas flow based on mass conservation required a shift in impact time of τ 0 = −5 ms to overlap with the bulk velocity obtained in both FOV's of the PIV measurements (Fig. 4). Additionally, a correction factor (α i ) is required to overlap the measured bulk velocity with the bulk velocity based on mass conservation. The correction factors are 1.54 and  5). Note that the velocity profile of FOV1 is different from that of FOV2 (Fig. 5a). First, the averaged valid velocity data of FOV1 is shown in 5 a. The measurements of FOV1 are not as repeatable as those of FOV2 with error bars that are smaller than the marker size (typically 9.5% of the velocity) for 0.3 ≤ x/∆x ≤ 0.75. However, the gas velocity decreases close to the impact wall and the wave tip. The gas velocity at the impact wall (x/∆x = 0) is lower than expected, as the spatial resolution of FOV1 is limited. The center of the domain of FOV1 is slightly misaligned with the center of the impact wall, which results in a small perspective error. The perspective error blocks the view of the tracer particles for approximately x/∆x ≤ 0.1. Furthermore, the detection of accurate velocity vectors in the neighborhood of walls is compromised and errors can easily propagate through the domain [59]. Additionally, loss of correlation occurs, which decreases the velocity at the impact wall for FOV1. The domain of FOV2 has been properly aligned with a higher spatial resolution, which does not result in a decrease of the velocity near the impact wall. The measurements of FOV2 are repeatable with error bars that are smaller than the marker size (typically 2.6% of the velocity) for 0 ≤ x/∆x ≤ 0.75.
On the other hand, the decrease in gas velocity at the wave tip (x/∆x = 1) is also observed for FOV2.
The vertical gas velocity is observed to decrease near the wave crest for approximately 0.75 ≤ x/∆x ≤ 1, which is denoted as the sub-window in panel a of Fig. 5. In panel b of Fig. 5 the same sub-window is shown for the smooth measurement (M226) of FOV1, which shows a non-uniform seeding distribution near the wave tip. The non-uniform seeding distribution consistently occurs over multiple measurements. The seeding is not only nonuniform, it also shows the displacement of streaks near the tip of the wave crest. The discontinuities in seeding concentration are already introduced at the nozzle of the seeding supply. The acceleration of the gas stretches patches of seeding to streaks, which results in a reduced PIV accuracy near the wave crest tip.
The discontinuities in image intensity that result from the non-uniform seeding result in a decrease of the vector quality [54]. A symmetric phase only filter is already applied to increase the peak detectability as described in Appendix B. However, the lack of identifiable seeding particles near the tip of the wave crest still results in a deterioration of the accuracy and a decrease of the magnitude of the velocity vectors in the gas. The measured bulk velocity of the gas at t − τ 0 = −15.0 ms is, without correction, approximately 4.3 and 4.5 m s −1 for FOV1 and FOV2 respectively. The bulk velocity for the gas based on mass conservation is approximately 5.2 m s −1 , which is significantly higher than the measured bulk velocity for both FOV's.
The velocity profile between the wave crest tip and the impact wall resembles that of a potential flow over a bluff body. Furthermore, the flow between the wave crest tip and the impact wall is in good approximation inviscid and irrotational, which enables a potential flow comparison. The Reynolds number based on the mass conservation bulk velocity at t − τ 0 = −15.0 ms is approximately Re = v b L/ν g ≈ 2 × 10 4 for a wall to tip distance of ∆x = 50 mm, and a kinematic viscosity of ν g = 15.06 × 10 −6 m 2 s −1 of air at 20 • C and standard atmospheric conditions. The measured velocity profile is fitted to the classical solution of a potential flow solution past a cylinder [60,61]. The proposed fit of the velocity profile is, based on the classical solution, presented as: v with the unknown dimensionless parameters C 0 , and C 1 , the velocity past the potential flow object (U 0 ), and the upwards velocity of the wave (v through ). In a previous manuscript the upwards velocity of the wave (i.e., the relative motion of the wave tip) was shown to be v through ≈ 1.23 m s −1 [43].
The potential flow object is fixed on the wave tip, which requires a translation from the laboratory coordinate system to the wave coordinate system. The trough velocity (v through ) is used to define the wave tip coordinate system. The remaining parameters of Eq.  The gas velocity near the wave crest is significantly higher with a velocity of 12 m s −1 compared to the bulk velocity of 5.2 m s −1 based on the fit. Recent numerical simulation also show the high gas velocity near the wave crest [62]. A higher gas velocity near the wave crest can result in an earlier onset of shear instability development (i.e., Kelvin-Helmholtz instability). The earlier onset of shear instabilities are, based on the fit, hypothesized to be the cause of the disturbances on the wave crest geometry, such as those shown in panel b of Fig. 3.

Dynamics of the circulation zone
The source of variability in peak impact pressure observed during wave impact experiments is still largely unknown. However, the formation of instabilities on the wave crest is hypothesized by others to be the main source of wave impact pressure variability [22,35]. Does the gas flow induce additional variability apart from the formation of instabilities? In this section the gas flow field over the wave crest is discussed. This includes a discussion on the development of large vortices over the wave crest and of local vortices near the wave crest tip.
The vorticity component perpendicular to the measurement domain (ω z = ∂v/∂x−∂u/∂y) is estimated with a filtered secondorder difference method [52]. This method is based on an 8-point estimate of the local circulation. The presence of the free surface locally complicates the calculation of the vorticity, as the vorticity cannot always be determined with an 8-point estimate [44]. A free surface that is included in the computation of the vorticity will result in an artificial shear-layer at or near the surface [63]. Therefore, the vorticity is only calculated with the valid vectors fully in the gas phase. In this section two waves will be highlighted that display differences in the development of their vorticity field, even though the wave shape and velocity are similar on a global scale. The two waves are representative of the larger set of measurements performed in this study. The first wave has been denoted as the perturbed measurement (M225), whereas the second wave will be denoted as the smooth measurement (M226). The perturbed measurement (M225) displays a continuous growth of the circulation zone over the wave crest, whereas the smooth measurement (M226) tends to separate prior to impact.
The vorticity development is shown in Fig. 6 for the perturbed wave (M225) a typical measurement with a disturbance. Initially for t = −24.0 ms, there is no visible vorticity development over the wave crest (Fig. 6a). If a vortex is present at this stage it cannot be resolved because of the limited vector spacing of 3.8 mm. However, a small vortex and shear-layer are visible at t = −18.0 ms (Fig. 6b). The vortex structure forms on the back of the wave crest close to the disturbance on the free surface [44,64]. Gas flow separation over free surfaces is complex to detect, but it is often linked to wave breaking during wind wave generation [64]. The vortical structure grows in size and the vorticity magnitude increases as the wave approaches the wall (Fig. 6). The vortex grows behind the disturbance, which becomes especially clear for t = −6.0 ms in panel d of Fig. 6.
A vortical structure consistently develops at the backside of the wave crest for all measurements. Two typical vortical structures can be identified in the current experiments based on the shape of the free surface, the growth rate and size of the vortex. The measurement with the perturbed wave (M225) is one of the typical examples of the first vortical structure, that separates from the wave crest close to the disturbance on the free surface. On the other hand, the measurement with the smooth wave (M226) is one of the typical examples of the second vortical structure. Fig. 7 shows the second typical vortical structure for the measurement with the smooth wave (M226). A clear difference between both waves and the development of the vortex is already present at t = −24.0 ms. The wave crest of the smooth wave (M226) is smoother than that of the smooth wave (M225) and this influences the vorticity development. For example, a small but consistent vortex is already present for t = −24.0 ms (Fig. 7a) that grows in both size and vorticity magnitude over time (Fig. 7b-c). Prior to impact, however, the vortex remains approximately stationary and equal in vorticity magnitude (t ≥ −12.0 ms). Additionally, a secondary vortex develops near the wave tip (Fig. 7d). The circulation (Γ ) of the vortex is further evaluated (Fig. 6d). The focus is especially on the typical behavior of the vorticity development for either a smooth wave crest (M226) or a wave crest that contains a disturbance (M225). The selected measurements highlight differences in the global vortex development for similar waves. The other waves of this set display similar variations in the global vortex development.  However, the differences in the global vortex development are more subtle.
A vortex identification technique is required to determine the geometric and vortical properties of the vortex (Fig. 6b). There are several vortex identification techniques available, many of which are based on derived quantities of the velocity gradient tensor [65]. Vortex identification techniques that are based on the velocity gradient tensor are not robust if a large vortex is superimposed on a turbulent velocity field [66]. On the other hand, the resolution of the velocity field determined over FOV1 is such that the small scale turbulence is not resolved, which inhibits the use of methods that rely on the flow field topology. The swirling strength criterion (λ ci ) is, therefore, used to identify the existence and center of a vortex [67]. The swirling strength criterion is frame independent and only identifies regions of circular motion. The swirling strength criterion is defined as the imaginary part of the complex eigenvalue pair derived from the eigendecomposition of the velocity gradient tensor [67].
The center of the vortex is defined by the maximum of the swirling strength criterion. The swirling strength criterion also defines the contours of the vortical structure. The contour of the vortical structure is defined by the values of the circulation strength that are higher than a specified threshold (i.e., λ ci > 2.5λ ci,rms ). Minor differences were observed for different thresholds, but no qualitative differences were identified. The rootmean square of the swirling strength criterion is only determined over the non-zero components. Thereafter, the circulation (Γ ) of the vortex is calculated as, where ω z is the vorticity, and N the number of elements inside the contour of the vortical structure defined by the swirling strength criterion.
The circulation (Γ ) is determined for all measurements of FOV1 (Fig. 8). However, only the data points for the perturbed The saturation of circulation in a vortex is a well-known phenomena [68]. An associated time scale (i.e., the formation number) can be defined as t * = Ut/L with U a typical velocity scale and L a typical length scale. The formation number is a dimensionless measure of the time it takes to create the largest possible vortex, which Gharib et al. [68] postulates occurs at approximately t * = 4. The circulation can be increased after the maximum formation number, but the vortex core will then shed the accumulated vorticity in its wake [68].
In Fig. 8 the maximum circulation is not obtained at a comparable time for the perturbed (M225) and smooth (M226) measurement, which is interesting as the global wave behavior for both waves is repeatable [43]. Typically a length and time scale are easily defined and relatively constant. For example, Gharib et al. [68] showed the formation of a vortex ring from the motion of a piston where the length and time scale of the flow are defined by the pistons shape and its motion. In the current work, the typical length and time scale are non trivial and probably not even constant over time. However, an estimate of the required length scale is obtained by assuming a formation number of t * = 4 at the maximum circulation of the perturbed measurement (M225). This can be obtained with a velocity over length scale ratio of approximately U/L ∼ 0.5. A typical velocity scale of the plunging breaking wave is the velocity of the wave tip at 3.31 m s −1 , which would imply an unphysical length scale of approximately 6.6 m. Consequently, the circulation data is not displayed against the formation number. However, the largest possible vortex remains interesting, as the vortex is also observed to separate from the wave crest. For example, the vortical structure remains approximately stationary for the perturbed measurement (M225) for t > −12.0 ms (Fig. 6c-d), which overlaps with the circulation limit (Fig. 8).
The different circulation limits between the perturbed (M225) and smooth (M226) measurement and the occurrence of secondary vortices in the smooth measurement (M226) warrant a detailed investigation of the velocity field over the wave crest tip. However, the resolution of the velocity field of FOV1 is not high enough for a detailed investigation of the circulation and vorticity near the wave crest tip. Therefore, the development of the wave crest tip and flow field are studied using the water droplets as a visualization tool (Fig. 9). Fig. 9 is also available as a movie (see Supplementary Material).
In panel a and b of Fig. 9 the visualizations of respectively the perturbed (M225) and smooth (M226) measurement are shown for t = −14.4 ms. The disturbance on the wave crest is directly visible for the perturbed (M225) measurement, whereas it is absent for the smooth measurement (M226). Furthermore, a part of the vortex wake is identified in panel a, whereas for the smooth measurement (M226) the vortex is further downstream and is not visible in the selected window. The smooth measurement (M226) displays the development of a small recirculation region near the wave crest tip.
The recirculation region near the wave crest tip develops further for both measurements at t = −12.0 ms (Fig. 9c-d).
A small secondary vortex is already visible for the perturbed measurement (M225), whereas it is still absent for the smooth measurement (M226). Interestingly, the secondary vortex develops further for the smooth measurement (M226) at t = −8.8 ms, whereas for the perturbed measurement (M225) the development of the secondary vortex is blocked by the upstream disturbance on the wave crest. The secondary vortex is observed to breakdown for the perturbed measurement (M225) at t = −8.8 ms and a shear layer develops over the wave crest. The breakdown of the secondary vortex in the perturbed measurement (M225) seems to correspond to the saturation limit of the circulation (Fig. 8). The secondary vortex remains attached to the wave crest tip for the smooth measurement (M226) (Fig. 9h). The development of a shear layer after separation of a vortex resembles that of a wing that experiences stall, which results in a decrease of the lift force and an increase in the drag force [69]. Fig. 10 shows the tip coordinate of respectively the perturbed (M225) and smooth (M226) measurement obtained from the stereo-PLIF measurements. The difference between both measurements is clear, especially, prior to impact where the y wtcoordinate of the wave tip displays a difference of approximately 6.8 mm. Initially, for t < −20 ms, both measurements show a similar displacement of the wave tip y wt -coordinate. The development of the recirculation zone near the wave crest tip results in an increase of the lift force for t > −12.0 ms in the smooth measurement (M226) (Fig. 9d). The lift induced by the secondary vortex result in an additional upward displacement of the thin wave tip.
Wave tip deflection is observed in numerous numerical simulations (see for example Behruzi et al. [58], Etienne et al. [62]). The simulations often initialize the global wave generation with a potential flow model. Thereafter, the computations are performed with both liquid and gas coupling, which -as we show -is essential to model the variability of wave impact pressure. Especially, the development of circulation is important when considering the impact pressure variability, as the vortex at the wave crest tip

Conclusion
The repeatability of the gas flow over plunging breaking waves is determined prior to impact. Several repeatable focused waves are generated that impact on a vertical wall. The global gas flow behavior is comparable to that of a plunging breaking wave [45] Fig. 10. The wave tip location for two typical measurements M225 with a disturbance on the wave crest and M226 with a smooth wave crest surface. The difference in wave crest shape results in a vortex with a lower circulation for the smooth (M226) measurement, whereas the sharp disturbance on the wave crest of the perturbed measurement (M225) results in a significantly higher recirculation zone. The formation of a secondary vortex near the tip of the wave for the smooth measurement (M226) results in a change in slope of the y wt coordinate. The circulation zone has a significant influence on the wave tip impact location, and in turn on the variability of the wave impact pressure. and a progressive wave [44]. A vortex develops over the leeward side of the plunging breaking wave that separates from the breaking wave and remains stationary in the stagnant air [44].
The gas flow drives the instability development on the wave crest, which requires an appropriate estimate of the bulk velocity. The bulk velocity was derived previously based on the shape of the global wave [43]. However, the integrated velocity profile between the wave tip and the impact wall results in a significantly lower bulk velocity compared to the mass conservation estimate. The high velocity close to the wave crest tip is not resolved, because of the limited resolution of the PIV measurements and the absence of sufficient seeding near the wave tip. A fit of the velocity data that resembles a potential flow estimate over a bluff body results in a significant improvement of the bulk velocity estimate. Consequently, the gas velocity close to the wave tip is approximately 2 times higher than the bulk velocity estimate, which may results in an earlier onset of instability development.
Two typical vortical structures are shown to develop over the wave crest for which the perturbed (M225) and smooth (M226) measurement are illustrative examples, that either have a smooth wave crest (M226) or a wave crest with a disturbance (M225). The development of the vortex is different for each typical measurement. The vortical structure becomes larger if a disturbance is present on the wave crest surface compared to the other typical measurement. The circulation magnitude is almost double that of a smooth wave crest. On the other hand, the circulation development is similar for both measurements, where the circulation remains stationary after an initial growth. The circulation development is typical for a vortex that eventually separates from the wave crest at a specified dimensionless formation number, but in this measurement a typical formation number cannot be defined [68].
The gas flow may affect the wave impact pressure variability, as it enhances the variability in the wave impact location. For example, a secondary vortex develops near the wave crest tip for the measurement of an initially smooth wave crest (M226). On the other hand, the disturbance on the wave crest for the perturbed measurement (M225) blocks the development of a secondary vortex, which results in a shear-layer with minimal lift. The secondary vortex increases the lift near the wave crest tip, which results in an additional deflection of the wave crest tip. The deflection of the wave crest tip is approximately 6.8 mm compared to the other typical measurement, which is significant on the scale of a typical pressure sensor. Consequently, an additional mechanism for the development of variability in the wave impact location and possibly impact pressure variability on globally similar waves is the development of a secondary vortex (i.e., flow separation) over the wave crest tip.
The global characteristics of an air-water wave impact are unaltered by the presence of the gas phase. The global characteristics of the wave impact pressure can be retrieved with pressure impulse models [15]. On the other hand, the variability of wave impact location and thereby probably the wave impact pressure is largely driven by local phenomena, such as flow separation and the development of free surface instabilities. Consequently, numerical models that aim to quantify wave impact pressure variability on a local scale require accurate models of the gas phase. Additionally, the variability in the wave impact location should also be studied at different atmospheric conditions (i.e., density ratios and ullage pressures) and different global wave shapes. A simultaneous study of the detailed wave impact dynamics (i.e., the wave tip location) and wave impact pressure should be conducted by others to confirm the variability in impact pressure induced by the variation in wave impact location and, thereby, confirm the importance of gas flow separation over the wave crest surface.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
The data that support the findings of this study are available from the corresponding author upon reasonable request.

Acknowledgments
This work is part of the public-private research program Sloshing of Liquefied Natural Gas (SLING) project P14-10. The support by the Netherlands Organization for Scientific Research (NWO) Domain Applied and Engineering Sciences, and project partners is gratefully acknowledged. The authors are grateful to the lab technicians at both the Laboratory for Aero-and Hydrodynamics and the Hydraulic Engineering Laboratory of the Delft University of Technology for their assistance with the experimental setup and their support during the experiments.  11. The probability density function (PDF) of the tracer particle size is shown. A typical log-normal fit is applied to the data with a mean of µ = 2.2 µm and a standard deviation of σ = 0.27 µm. The inset shows a typical particle image with a normalized interference pattern. Fig. 12. The correlation peak is determined at a typical location near the wave crest tip for FOV1 at the final pass of the PIV algorithm. (a) The correlation peak is typically determined with a cross correlation PIV algorithm. (b) The correlation peak is determined with a symmetric phase-only filter. The correlation peak is clearly defined when the image data is filtered with a symmetric phase-only filter, whereas it is obscured by the particle streak in panel a.

Appendix A. Interferometric particle imaging
An interferometric particle imaging (IPI) technique is used to determine the tracer particle size [70][71][72]. The particles are generated with a commercial ultrasonic humidifier (UHW, Medisana). The tracer particles are illuminated with a 1 mm thick laser light sheet generated by a high-speed Nd:YLF laser (LDY 304 PIV laser, Litron) with a perpendicular polarization and a wavelength of λ = 532 nm. Images are acquired with an high-speed CMOS camera (Fastcam APX-RS, Photron) equipped with a 200 mm Micro-Nikkor objective at an f-number of f # = 4 that results in an aperture for the lens system of d a = f /f # = 50 mm. The camera is initially focused on a two-level double sided calibration plate (Type 58, LaVision) to obtain a magnification of M 0 = 1.2 at a distance of z = 424 mm. After calibration the optical system is slightly defocused to obtain interference images of the particles.
The camera is placed at an angle of θ r = 90 • with respect to the laser light sheet to increase the signal to noise ratio for the selected polarization [70]. There is a direct relation between the number of fringes (N) and the particle size (d p ). For a refractive index ratio (m = n l /n g > 1) and an aperture size that is small compared to the focal distance (d a /z ≪ 1) the relation is [71]: [ cos(θ r /2) + m sin(θ r /2) √ m 2 − 2m cos(θ r /2) with m = n l /n g = 1.33/1 the ratio of the refractive index of water and air. The particle size is solely determined by the optical choices in the set-up, such as the focal distance and the aperture size. The viewing angle only defines the signal to noise ratio of the measured particle images. However, the IPI technique is not well suited for dense sprays. In the current work, the spray was diluted to obtain usable particle images. The PDF of the particle size is shown in Fig. 11, where the distribution has an average particle size of d p = e µ+σ 2 /2 ≈ 9.8 µm.

Appendix B. Symmetric phase only filtering
The symmetric phase only filter is used, with a window function defined by Wernet [54]: with |S i | the real part of the Fourier transform. The Fourier transform of the interrogation window located at coordinates p, q is defined as S i (p, q) = |S i (p, q)|e −iφ(p,q) for camera i. The window function (i.e., filter) is only applied to a single image Fourier transformed of the image pair. The correlation peak of the symmetric phase only filter is compared to a typical correlation method in Fig. 12. A single interrogation window is selected for a typical PIV image. The correlation peak is determined with both a normal cross correlation (Fig. 12a) and a phase only filtered correlation (Fig. 12b). The correlation peak is clearly identified for the phase only filtered window, whereas it is unclear in the original image.