In this work, a detailed description of the internal flow field in a collapsing bore generated on a slope in a wave flume is given. It is found that in the case at hand, just prior to breaking, the shape of the free surface and the flow field below are dominated by capillary effects. While numerical approximations are able to predict the development of the free surface as it shoals on the laboratory beach, the internal flow field is poorly predicted by standard numerical models.

## I. INTRODUCTION

Wave breaking is an important and ubiquitous phenomenon which happens virtually in all flows involving a free surface. As clearly brought out by a number of review articles and monographs,^{3,6,18,29,31,38,39,42,43} wave breaking has been a major focus of wave and ocean research for a long time. While it is well understood that wave breaking is central for the study of the energy budget of the oceans and air-sea interaction, it is also noted in the above works that the nature of wave breaking remains poorly understood. Indeed, wave breaking is a classical multiple-scale problem which exhibits a number of complicating factors in the flow such as circulation, turbulence, capillarity, and intermittency.

Wave breaking happens in many shapes and forms, and different types of breakings require a variety of methods for study. The object of the present paper is the study of a collapsing breaker on a moderate slope. Collapsing breakers appear on the margin between plunging and surging breakers, as found in the work by Galvin.^{22} Possibly because of their relatively rare occurrence, studies of collapsing breakers are few and far between. Indeed, collapsing breakers are not explicitly mentioned in Ref. 25, where the authors define a detailed classification of breaking waves on a planar beach in terms of a single parameter *S*_{0} based on the bottom slope, waveheight, and wavelength.

The focus of the present work is twofold. On the one hand, using new experiments performed in a wave flume at National Chung Hsing University in Taiwan, we aim to give a detailed description of the internal flow including velocity, Lagrangian acceleration, and pressure fields in the prebreaking stage of a capillarity-dominated collapsing breaker. In particular, we aim to describe details of the fluid flow below the free surface in order to identify a number of indicators for the onset of wave breaking.

Second, the internal flow field is compared with the results of numerical simulations using a Boussinesq model.^{48} It is found that this numerical model has certain challenges when it comes to describing the internal flow field in a flow of this kind.

In order to create a collapsing breaker, a solitary wave of a certain waveheight is initiated with a wavemaker at one end of the wave flume. The wave propagates through the tank and starts deforming as it comes upon the slope located at the point *x* = 0. As the wave steepens, the wavefront starts to resemble a bore rushing onto the slope and eventually collapses without spilling or overturning. In this laboratory scale experiment, in the final stage before the bore collapses, the flow near the free surface is dominated by capillary effects, which leads to interesting pressure anomalies. The shape of the bore just before breaking takes a *dolphin head*-like appearance, resembling similar profiles found in the context of capillary jumps such as discussed in Refs. 20 and 37. The shape of the wavefront is also similar to a spilling breaker in the presence of strong capillarity, such as discussed in Refs. 18 and 19.

The collapsing of the breaker is eventually indicated by the development of negative vertical velocity components which are necessitated by the acceleration due to the pressure anomalies resulting from strong capillarity. As the bore collapses, its internal flow structure resembles a plunging breaker (see, for example, Ref. 49) even though the free surface is still intact in the early stages of wave breaking. Eventually, as the toe of the capillary region gets ever more acute, it acts essentially as a hydrofoil, spoiling the free surface and creating a vortical motion in the breaking bore. The dynamics of the ensuing eddy motion could possibly be described numerically, but this is beyond the scope of the present article. Indeed, it is not clear whether existing point-vortex models such as, for example, those used in Ref. 15 will be able to capture the surface and eddy motion to a satisfactory degree of accuracy.

## II. EXPERIMENTAL TECHNIQUES

### A. Experimental apparatus

The experiments were conducted in a wave flume located at the Department of Civil Engineering, National Chung Hsing University, Taiwan. The internal dimensions of the wave flume are 14.00 m long, 0.25 m wide, and 0.50 m deep. The bottom and two sidewalls of the wave flume were made of tempered glass to allow optical access. A piston-type wavemaker driven by a servo-motor is mounted at one end of the wave flume. The method proposed by Goring in Ref. 23 for generating the solitary wave is used. As shown in Refs. 33–36, highly repeatable solitary waves can be produced by this wavemaker.

A sloping beach made of tempered glass with a slope of 1:20 was installed in the wave flume. This is a moderately steep slope, much more gentle than the steep slopes used in Refs. 28 and 52. The toe of the sloping beach was fixed 6.48 m away from the wave paddle at rest. A Cartesian coordinate system with the origin (*x*, *y*) = (0, 0) cm being located at the toe of the sloping beach is used. The x-axis is oriented in the horizontal direction and measured positive onshore from the toe. The y-axis is oriented in the vertical direction and measured positive upward from the horizontal bottom. The schematic diagram of the experimental setup and coordinate system is shown in Fig. 1.

The free surface elevation was measured using two ultrasonic wave gauges. The first gauge located at *x* = −150 cm was used to measure the time series of free surface elevation *η*_{0} and the waveheight *H*_{0} for the incident solitary wave propagating over the horizontal bottom with a still water depth of *h*_{0}. A voltage signal output from this wave gauge was employed to trigger the camera for capturing the images. Furthermore, the second gauge is placed at the toe at *x* = 0 cm to precisely identify the time, *t* = 0 s, at which the wave crest exactly reaches the toe. As shown in Fig. 1, *η*(*x*, *t*) is the instantaneous wave profile and *h*(*x*) is the still-water depth at a specified location over the sloping beach.

### B. Flow visualization and velocity measurements

The structure of the prebreaking, breaking, and postbreaking waves was explored using flow visualization techniques (FVT) and high-speed particle image velocimetry (HSPIV) measurements. Titanium dioxide (TiO_{2}) was used for the seeding particles. These seeding particles have a refractive index of 2.6 and a mean diameter of 1.8 *μ*m, together with mean specific gravity and settling velocity (estimated by Stoke’s law) equal to 3.547 and 4.5 × 10^{−4} cm/s, respectively. The fall velocity is very small and thus ignored as compared with the typical velocity of interest in this study. An argon-ion laser head (Innova-300, Coherent Inc.) with a maximum power of 4 W was used as a light source. A fan-shaped laser light sheet (1.5 mm thick) was formed while the laser beam was guided by three reflection mirrors and then passed through the cylindrical lens. The light sheet, located 8 cm away from the glass sidewall, was eventually projected upward through the glass bottom of the wave flume to illuminate the seeding particles suspended uniformly in the water column. A high-speed digital camera (Phantom M310, Vision Research) with a maximum framing rate of 3260 Hz under the largest resolution of 1280 × 800 pixels was used to capture the images of both the free surface profile and the flow structure underneath the free surface.

The flow visualization images were captured by using a particle-trajectory photography method, with a high exposure time to allow the path-line of the particles to be captured. The framing rate and the exposure time of the high-speed camera for taking flow visualization images were 500 Hz and 1900 *μ*s. On the other hand, to ensure high time-resolved HSPIV measurements, the framing rate was 2000 Hz, much higher than FVT, and the exposure time was 490 *μ*s, smaller than FVT. The images of flow fields with the illuminated moving seeding particles were recorded continuously by using the high-speed camera. The free surface profile of the shoaling solitary waves, *η*(*x*, *t*), was extracted from the images of the HSPIV measurements. The time interval in between each such image was very short, and a high-resolution record of *η*(*x*, *t*) could be obtained. Furthermore, due to the mirror effect of the interface between the air and water, the particles on the free surface were illuminated by the laser light sheet and can be clearly identified.

One field of view (FOV1) was employed in this study. As shown in Fig. 1, the FOV1 was located between *x* = 270.76 cm and *x* = 276.34 cm, *y* = 13.18 cm and *y* = 16.66 cm. The width and height of FOV1 were 5.58 cm and 3.48 cm, respectively. The resolution was 4.356 × 10^{−3} cm/pixel. A cross-correlation calculation was used to determine the velocity vector from the images in which the brightness of the seeding particles had been intensified by using the Laplacian edge-enhancement technique.^{1} The multipass PIV algorithm was then employed to calculate the instantaneous velocity field from three correlated images. The calculation is started from an interrogation window size of 64 × 64 pixels and ending with a window size of 8 × 8 pixels with a 50% overlap. Consequently, spurious vectors were removed by employing both global-range and median filters. Missing vectors are then interpolated to complete the whole velocity vector field.

### C. Experimental conditions

An incident solitary wave with waveheight *H*_{0} = 1.12 cm was created on water of an undisturbed depth of *h*_{0} = 14 cm so that the ratio $H0\u2032=H0/h0$ was equal to 0.08. The wave conditions are listed in Table I. Ten repeated runs were conducted for the HSPIV measurements. An ensemble average of all runs was used to describe the spatial or temporal variation of the velocity fields.

H_{0} (cm)
. | h_{0} (cm)
. | H_{0}/h_{0}
. |
---|---|---|

1.12 | 14.0 | 0.080 |

H_{0} (cm)
. | h_{0} (cm)
. | H_{0}/h_{0}
. |
---|---|---|

1.12 | 14.0 | 0.080 |

In order to obtain a collapsing breaker, the classification of Grilli *et al.*^{25} was used. In that work, the authors define the parameter

where *s* = tan *θ* is the slope, *h*_{0} is the undisturbed depth prior to the sloping bottom, and *L*_{0} is the wavelength. Following Ref. 25, the wavelength in Boussinesq’s solitary wave theory is measured at the point of maximal slope on the wave profile, which leads to the relation

With this definition of *L*_{0}, *S*_{0} is given in terms of the slope *s* and relative waveheight $H0\u2032$ by

Grilli *et al.*^{25} observed the following breaker types:

spilling:

*S*_{0}< 0.025,plunging: 0.025 <

*S*_{0}< 0.3, andsurging: 0.3 <

*S*_{0}< 0.37.

In the present case, we are using a slope of 1:20 so that *s* = 0.05. Since the relative waveheight is $H0\u2032=0.08$, the parameter value for the current experiment is *S*_{0} = 0.269. Thus, the waveform used here is in the range of plunging breakers, but close to surging according to the classification above. As noted by Galvin,^{22} some waves fall between plunging and surging and can be classified as collapsing breakers, and this is the case here.

## III. DESCRIPTION OF THE FLOW

### A. Wave generation and propagation

A solitary wave was generated in the wave flume by executing a quarter stroke with the piston wavemaker. The early stage of the evolution of the solitary wave was checked against a potential flow solver,^{24,26} where the solitary wave was generated numerically using a variant of Tanaka’s method (see Refs. 14 and 51). The resulting time series can be seen in Fig. 2. Note that in both time series, there is a nearly perfect match between experimental data and numerical simulation.

### B. Details of the wave breaking process

As the wave passes the toe of the sloping bottom, it starts to feel the upward slope and the wave starts to steepen. The slope provides a force in the direction opposing the wave motion, which eventually stops the motion. The reflection then leads to a rundown and backward motion. The rundown is not in focus in the present study.

As the wave collapses on the beach, there is little air entrainment. Moreover, the width of the tank is relatively small at 25 cm. With these restrictions, we notice that the flow is not being able to form three-dimensional turbulent structures even in the postbreaking stage. In a sense, the flow is what one might call quasilaminar.

Next, we describe the details of the prebreaking development. From Fig. 3, it can be seen that fluid velocities near the top of the surging wave become larger. The wavefront gets steeper as more and more fluid concentrates behind the surging wavefront, creating a bulging “dolphin head” shape. A study of the individual terms in the Lagrangian acceleration, i.e., the terms $\u2202u\u2202t$, $u\u2202u\u2202x$, and $v\u2202u\u2202y$, shows that the large velocities in Fig. 3 are not due to acceleration of fluid particles, but rather due to convective effects and concentration of fast particles behind the wavefront. In fact, consulting Fig. 4 shows that there is a fairly strong deceleration of particles behind the wavefront. As will be explained later, this deceleration is due to a large excess pressure just behind the wavefront due to capillary effects. The same capillary effects also contribute to acceleration of fluid particles in the free surface (cf. Fig. 5).

The onset of breaking may be defined as the first time the vertical velocity component of a particle behind the wavefront is negative. In Fig. 7, it can be seen that the vertical velocity component of particles in a larger and larger region starts to turn negative. This incipient downward motion constitutes the beginning of the creation of an internal circulation behind the wavefront.

In Fig. 5, particles in the free surface are being followed (manually), and it can be seen that these particles are accelerating while the free surface is decelerating. This feature seems to be common in waves approaching breaking. For example, similar behavior was found recently in deep-water waves which approach the breaking point.^{30} It is also interesting to view these data in light of recent work^{7} where a universal breaking criterion was put forward. This criterion can be formulated in terms of the horizontal component of the fluid particle velocity at the crest *u* and the velocity of the crest itself *c*, and it states that a wavetrain is liable to feature wave breaking if the ratio *B* = *u*/*c* exceeds the threshold 0.85–0.86, which is in contrast to the usual convective criterion (see Ref. 12 and references therein) which places the critical value at about 1. In the present case, as the wave enters the slope, it forms a steep front so that there is no well defined wave crest. Nevertheless, using the data shown in Fig. 5, we may compare the front velocity to the horizontal component of the particle velocity. A representative result is shown in Table II. Note that the value *B* ∼ 0.87 is achieved about 0.003 s before the wave starts to break (cf. Fig. 7).

t (s) | 2.656 | 2.661 | 2.666 | 2.671 | 2.676 |

B | 0.75 | 0.77 | 0.78 | 0.84 | 0.87 |

t (s) | 2.656 | 2.661 | 2.666 | 2.671 | 2.676 |

B | 0.75 | 0.77 | 0.78 | 0.84 | 0.87 |

Finally, let us explain the role of capillarity in the prebreaking development of the wave. The horizontal momentum balance is written in terms of the horizontal component of the velocity field *u* and the stress vector *σ*_{x} as

where $DDt$ represents the material (Lagrangian) derivative and *ρ* is the fluid density. Given that the fluid velocity near the wavefront is similar to the velocity of the wave itself, we do not expect boundary layer effects to be dominant at the free surface so that the momentum balance reduces to

where *p* is the fluid pressure. Now from Fig. 4 which shows the Lagrangian acceleration in the *x* direction, we see that $DuDt$ is negative in the upper part of the wave. Thus, according to the above formula, the pressure *p* increases in the direction of wave propagation in the upper part of the wave. In contrast, in the lower part of the wave (red part), the Lagrangian acceleration $DuDt$ is positive so that the pressure decreases in the direction of wave propagation.

These findings can be explained by looking at the free surface condition with capillarity. Indeed, the balance of forces is written in terms of the fluid pressure *p*, the atmospheric pressure *p*_{a}, the free surface excursion *η*, and the capillary parameter *τ* as

This formula clearly shows that near the head of the wave, where the free surface is convex, the fluid pressure is above atmospheric (i.e., the gauge pressure *p* − *p*_{a} is positive), while near the bottom of the wave, the fluid pressure is below atmospheric since the free surface is concave.

Thus, as indicated in Fig. 6, there is an acceleration of fluid particles on the free surface, while particles behind the leading front of the wave are decelerated. Moreover, together, these pressure conditions also lead to acceleration of fluid particles near the middle of the leading front in the negative *y*-direction. After carefully analyzing the data using HSPIV method, it was indeed found that the internal velocity field under the head of the wave develops negative vertical velocity components (see Fig. 7).

## IV. COMPARISON WITH THE BOUSSINESQ MODEL

As there are no closed-form solutions of solitary waves propagating and breaking on a slope, the experimental data will be compared to simulations done with a Boussinesq model. There are some theoretical works providing nearly closed-form solutions for standing waves on a slope (see Refs. 10, 13, and 50), but in the present situation, numerical simulation is the only reasonable choice.

### A. Model description

The Boussinesq Ocean and Surf Zone (BOSZ) model is a phase-resolving Boussinesq-type model for the computation of nearshore waves, wave-driven currents, infragravity oscillations, and ship wake waves (see, for example, Refs. 16, 32, 46, and 47). The governing equations are based on a conserved variable formulation of Ref. 41. The numerical solution handles the nonlinear shallow water part of the governing equations with a finite-volume scheme based on a total-variation diminishing (TVD) reconstruction method of up to 5th order and a Riemann solver. This combination ensures robust and accurate computation of fast flows over irregular terrain including moving boundaries (wet/dry cell interfaces). The frequency dispersion terms are based on a central-difference scheme of second order. Time integration is carried out with an adaptive Runge-Kutta time-stepping scheme, allowing up to 4th order accuracy. For most computations, such as this test case, a 2nd order time integration is sufficient; however, some problems with more dispersive waves require at least a 3rd order integration scheme. Due to the presence of space-time derivatives of the evolution variables in the momentum equations, systems of equations have to be solved to extract the flow speed at the end of each time step. The two systems are directionally independent of each other with data-dependencies arising only in the *x*- or the *y*-directions, respectively. The bottom friction is accounted for through the widely used Manning-Strickler formula based on a roughness coefficient, which represents the surface property of the experimental layout. Here, we choose *n* = 0.012*s*/*m*^{1/3} to match the smooth laboratory slope.

The input waves are generated at the left boundary in various forms. For the present case, the solitary wave was generated in two ways: (a) with an analytical solution at the boundary and (b) by feeding a time series from the wave gauge through the boundary. The flow velocity at the boundary is set by long wave theory in the form $u=\eta g/h0$. As shown in Fig. 8, both methods have led to near-perfect replication of the input wave conditions.

### B. Reduction of dispersion based on the free-surface Froude number

As the flow depth becomes very shallow on the slope, wave breaking has to be incorporated locally into the Boussinesq model. As Boussinesq models do not have an inherent wave-breaking mechanism, a numerical criterion needs to be used in order to maintain stability. The strategy used here is based on the *free surface Froude-number Fr*, which can be determined from the free surface flow velocities. The governing equations allow for calculation of the flow velocity at any position in the water column based on the horizontal velocities and under the assumption of a predescribed quadratic velocity profile used in the derivation of the Boussinesq system.

We know that for *Fr* > 1, the flow becomes supercritical and bores can develop. Therefore, if the free surface Froude-number exceeds 1.0, the dispersion terms are locally and momentarily set to zero, i.e., the solution reduces to the hyperbolic nonlinear shallow-water equations in the cells around the wave breaking front. The local *Froude number* at the free surface is determined by

where the velocity at the free surface is given by the quadratic profile, which arises from the Taylor series expansion of the horizontal velocity which is usually employed in Boussinesq type models.^{41} Since the model code checks for the *Fr*-number criterion at each time step, the solution at the wavefront can dynamically switch between the full Boussinesq-type solution and the hydrostatic nonlinear shallow-water equations. Similar methods for the detection and treatment of wave breaking were tested in Refs. 4, 5, 9, and 27 and compared to data from Refs. 21, 50, and 55. In particular, this is in contrast to the typical wave steepening observed in hyperbolic models.^{44} One may also use the quadratic ansatz used in many Boussinesq models to reconstruct the velocity profile at any point in the fluid column (see Ref. 11, for example), though better results would be expected using higher-order reconstructions (see Ref. 54) or full Euler or Navier-Stokes equations.^{26,8}

### C. Comparison

Considering the results over the small PIV window, the BOSZ model computes the shape of the collapsing bore reasonably well. Particular snapshots show very good agreement (see Fig. 9). The wave steepness depends mostly on the grid spacing. Δ*x* = 0.50 cm shows good agreement, whereas Δ*x* = 0.75 cm and Δ*x* = 0.25 cm lead to a slightly gentler and steeper bore front. In more detail, note that within FOV1 shown in Fig. 9, the BOSZ code essentially solves a shallow-water system due to the breaking criterion switching off the dispersive parts. While the numerical solver does not feature molecular viscosity, it features numerical dissipation, and it appears that with a grid size of Δ*x* = 0.75 cm, there is just the right amount of dissipation to slow the hyperbolic steepening so that very good agreement with the experimental data is obtained (see the leftmost curve in the lower panel of Fig. 9). On the other hand, the numerical model does not incorporate capillarity, which is the reason for the poor agreement with the experimental data in the rightmost curve in the lower panel of Fig. 9. With decreasing grid size, Δ*x* = 0.5 cm and Δ*x* = 0.25 cm, the numerical dissipation weakens and more pronounced hyperbolic steepening is observed in the numerical solutions. This choice yields a simulation which is closer to the hyperbolic nature of the equations to be solved, but the comparison with the experimental data is not as good (see the leftmost curves in upper and middle panels of Fig. 9). However, this choice leads to a better comparison with the capillary region further up the slope as can be seen in the rightmost curves in the upper and middle panels of Fig. 9. The velocity field can be reconstructed using the quadratic ansatz inherent in the Boussinesq model. An example is shown in Fig. 10. As can be seen in this Fig. 10, the intricate features of the velocity field due to capillarity at the bore front are not captured by the Boussinesq model. This is to be expected as capillarity is not part of the Boussinesq model since it is generally a long-wave model.

### D. Wave runup

The measured maximum run-up height for 14 cm still water depth and waveheight to water depth ratio being 0.08 is *R* = 3.154 cm. With *n* = 0.012*s*/*m*^{1/3}, the computed maximum runup agrees exactly with the measured runup of 3.154 cm. The computed maximum runup is the highest cell on the straight slope that was wet at any time during the computation. This corresponds to a minimum water depth of 0.001 cm.

## V. CONCLUSION

In a certain sense, one may view collapsing breakers as a hybrid between a plunging and a surging breaker. In a surging breaker, the steep slope provides a strong force in the negative *x*-direction (i.e., in the direction opposed to the propagation of the wave). This force is indeed strong enough to slow the wave, eventually arresting the wave motion without breaking.

In a spilling breaker, the slope is comparatively gentler, the force in the negative *x*-direction comparatively smaller, and the waveheight comparatively larger. As a result, the fluid motion in the lower half of the wave is inhibited while the force distributed through the fluid from the sloping bottom is not strong enough to arrest the fluid motion in the upper half of the wave. Thus, the upper half of the wave overtakes the lower half, leading to the well known formation of a jet. As the leading part of the jet feels gravity and starts to fall, a reconnection with the lower part of the wave is formed, and wave breaking ensues.

The collapsing breaker appears to represent a balance where the pressure forces provided by the slope are strong enough to slow the fluid motion in the whole of the wave, but not strong enough to prevent it from breaking. As a result, one may observe something which resembles an “internal plunging breaker.” This observation is most apparent when looking at camera footage from the PIV system, but can also be observed in Fig. 7.

Since a delicate balance between the bottom slope and the waveheight is required, collapsing breakers are not widely observed in the field. Nevertheless, there are accounts of collapsing breakers in the literature. For example, a photograph of a collapsing breaker is given on page 56 in Ref. 45. The importance of capillarity is usually gauged by looking at the Bond number *ρgL*^{2}/*τ*, for a length scale *L*. For large breakers such as the one shown in Ref. 45, the Bond number is large, and capillarity is not expected to play a decisive role.

On the other hand, the distinctive features of a collapsing breaker at the laboratory scale such as those discussed in the present paper are provided by capillarity. As shown in Fig. 4, the Lagrangian acceleration in the *x* direction is negative near the head of the wave and positive in the lower part. Using the horizontal momentum balance, this finding can be explained by capillary effects on the fluid pressure. In particular, capillarity contributes to excess pressure at the leading edge of the wave and a pressure deficit at the wave toe. Such pressure anomalies can sometimes be achieved in the presence of strong shear (cf. Ref. 2), but in the present case, capillarity appears to be the dominant effect.

The shape of the free surface of the collapsing breaker can be reasonably well approximated with a Boussinesq model. The BOSZ model used in this study was able to predict the main features of the free surface profile. On the other hand, the omission of capillary effects and the reliance of a quadratic velocity profile limit the applicability of the Boussinesq model to simulate the internal flow structure of the collapsing breaker at a laboratory scale.

From an operational point of view, capillarity has not been shown to be a major factor in coastal dynamics or the development of beach morphology, but it may be important in small-scale wave breaking, for example, in small parasitic breakers riding on the top of larger waves (see Ref. 40 and the references contained therein) or in smaller breakers occurring in the uprush on the beach.

The wave breaking in the experiment reported in here has been shown to occur at values of the breaking onset parameter *B* as defined in Ref. 7 which are close to the values indicated in Refs. 7 and 17. As the breaking criterion defined in these papers has been tested in deep and intermediate depths, it will be interesting to conduct further work exploring the relation between the breaking onset parameter *B* and the details of the wave breaking process in a larger number of cases in shallow water, and, in particular, in flows dominated by capillarity.

It should be noted that capillarity also comes to the fore in breaking waves at small scales in various other settings. For example, wave breaking in a two-phase flow in circular pipes has been studied recently in Ref. 53. Another interesting further problem will be to understand the generation of vorticity in the wave breaking process and the ensuing eddy motion.

## ACKNOWLEDGMENTS

This research was supported in part by the Research Council of Norway and the European Union’s Horizon 2020 research and innovation programme. P.G. was supported in part by the National Science Foundation, USA (Grant No. DMS-1615480). V.R. acknowledges financial support from the Isite program Energy and Environment Solutions (E2S), the Communauté d’Agglomération Pays Basque (CAPB), and the Communauté Région Nouvelle Aquitaine (CRNA) for the chair position HPC-Waves. J.-M.Y. was supported in part by the Ministry of Science and Technology, Taiwan (Grant Nos. MOST 106-2115-M-126-003 and MOST 107-2115-M-126-004). W.-Y.W., C.L., and M.-J.K. were supported in part by the Ministry of Science and Technology, Taiwan (Nos. 105-2221-E-005-033-MY3 and 106-2221-E-005-045-MY3).