W M O S T E R T

S P O P I N

L D E I K E

High-resolution direct simulation of deep water breaking waves: transition to turbulence, bubbles and droplets production

, with a distribution shape close to N d (r d ) ∝ r -2 d . The evolution of the droplet statistics appears controlled by the details of the impact process and subsequent splash-up. We discuss velocity distributions for the droplets, finding ejection velocities up to four times the phase speed of the wave, which are produced during the most intense splashing events of the breaking process.

Introduction

The broader context

The action of breaking waves on the ocean surface has a large and incompletely understood effect on the dynamics of mass, momentum and energy transfer between the ocean and the atmosphere, converting much of the wave energy into heat in a complex process which spans a wide range of scales [START_REF] Melville | The role of surface-wave breaking in air-sea interaction[END_REF]. Breaking also marks a transition at the ocean surface from laminar flow to two-phase turbulent mixing at small scales, modulating the dynamics of the upper ocean sub-mesoscales, particularly via Langmuir turbulence and fronts [START_REF] Mcwilliams | Submesoscale currents in the ocean[END_REF], and affects the transport of particles with implications for the fate of oil spills and plastic pollutants [START_REF] Deike | Lagrangian transport by breaking surface waves[END_REF]Pizzo et al. 2019). Furthermore, surface breaking injects a large amount of gas into the ocean via the entrainment of bubbles, including approximately 30% of the CO 2 which has been released into the atmosphere [START_REF] Deike | Gas transfer by breaking waves[END_REF][START_REF] Reichl | Contribution of sea-state dependent bubbles to air-sea carbon dioxide fluxes[END_REF]; breaking also ejects spray into the atmosphere, where it can convect and evaporate to leave salt crystals which may serve as cloud condensation nuclei [START_REF] De Leeuw | Production flux of sea spray aerosol[END_REF][START_REF] Veron | Ocean spray[END_REF].

Wave breaking involves transition from two-dimensional (2D) laminar wave flow to three-dimensional (3D) turbulence. As wave energy focuses through linear or nonlinear processes, local conditions on a wave surface become unstable and cause breaking, which transfers energy and momentum to the water column. The geometry and kinematics of the breaking waves have been extensively studied [START_REF] Longuet-Higgins | The deformation of steep surface waves on water-i. a numerical method of computation[END_REF][START_REF] Perlin | Breaking waves in deep and intermediate waters[END_REF][START_REF] Schwendeman | Sharp-crested breaking surface waves observed from a ship-based stereo video system[END_REF][START_REF] Fedele | Crest speeds of unsteady surface water waves[END_REF], and the identification of a breaking threshold with approaches based on the wave kinematics, dynamics or geometry remains a longstanding issue [START_REF] Melville | The instability and breaking of deep-water waves[END_REF][START_REF] Banner | Wave breaking onset and strength for two-dimensional deep-water wave groups[END_REF][START_REF] Perlin | Breaking waves in deep and intermediate waters[END_REF] with recent work discussing the link between the breaker kinematics and dynamics (Saket et al. 2017a;[START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF][START_REF] Pizzo | Theory of deep-water surface gravity waves derived from a lagrangian[END_REF].

While the initiation of the breaking phenomenon and the turbulence generated by it have been characterized [START_REF] Rapp | Laboratory measurements of deepwater breaking waves[END_REF][START_REF] Duncan | Gentle spilling breakers : crest profile evolution[END_REF][START_REF] Tulin | Laboratory observations of wave group evolution, including breaking effects[END_REF][START_REF] Melville | The velocity field under breaking waves: coherent structure and turbulence[END_REF][START_REF] Banner | Wave breaking onset and strength for two-dimensional deep-water wave groups[END_REF][START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF][START_REF] Drazen | Turbulence and mixing in unsteady breaking surface waves[END_REF], the time and length scales of the transition process remain to be explored. During this transition to turbulence, air is entrained, and bubbles are formed [START_REF] Lamarre | Air entrainment and dissipation in breaking waves[END_REF][START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF] and spray droplets are ejected [START_REF] Erinin | Spray generation by a plunging breaker[END_REF]. The measurements of 3D two-phase turbulence in the laboratory and in the field present many technical challenges in terms of successfully accessing the turbulent flow field and the size distributions of drops and bubbles during the active time of breaking.

Direct numerical simulations (DNS) therefore appear as an appealing tool. Owing to the computational difficulty and expense of modelling 3D multiphase flows, numerical studies began by using 2D breakers as analogues for the full 3D processes [START_REF] Chen | Two-dimensional navier-stokes simulation of breaking waves[END_REF][START_REF] Song | A numerical study of breaking waves[END_REF][START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF][START_REF] Iafrati | Energy dissipation mechanisms in wave breaking processes: spilling and highly aerated plunging breaking events[END_REF][START_REF] Deike | Capillary effects on wave breaking[END_REF]. Early development of nonlinear potential flow models has shed light on the breaking process up to the moment of impact [START_REF] Longuet-Higgins | The deformation of steep surface waves on water-i. a numerical method of computation[END_REF][START_REF] Dommermuth | Deep-water plunging breakers: a comparison between potential theory and experiments[END_REF], while 3D simulations have used reduced models such as large-eddy simulation to capture the breaking process itself [START_REF] Watanabe | Three-dimensional vortex structures under breaking waves[END_REF][START_REF] Lubin | Numerical simulations of three-dimensional plunging breaking waves: generation and evolution of aerated vortex filaments[END_REF][START_REF] Hao | Wind-wave coupling study using les of wind and phaseresolved simulation of nonlinear waves[END_REF], but the complete resolution of the breaker in DNS in 3D has only recently become feasible [START_REF] Fuster | Numerical simulation of droplets, bubbles and waves: state of the art[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF][START_REF] Wang | High-fidelity simulations of bubble, droplet and spray formation in breaking waves[END_REF][START_REF] Yang | Direct numerical simulation of wind turbulence over breaking waves[END_REF]. Surprisingly, despite the essentially 3D nature of the turbulence resulting from breaking, 2D breakers at the tested conditions have provided a reasonable estimate of the dissipation rates obtained from experiments and 3D computation (discussed further below). In contrast, the turbulent dissipation in internal wave breaking has been shown to be a clear 3D process [START_REF] Gayen | Turbulence during the generation of internal tide on a critical slope[END_REF].

Laboratory experiments and direct numerical simulations of canonical breaking waves

Canonical breaking waves have been studied using a variety of different approaches both experimental and numerical [START_REF] Duncan | An experimental investigation of breaking waves produced by a towed hydrofoil[END_REF][START_REF] Melville | The instability and breaking of deep-water waves[END_REF][START_REF] Rapp | Laboratory measurements of deepwater breaking waves[END_REF][START_REF] Melville | Energy dissipation by breaking waves[END_REF][START_REF] Duncan | Gentle spilling breakers : crest profile evolution[END_REF][START_REF] Banner | Wave breaking onset and strength for two-dimensional deep-water wave groups[END_REF][START_REF] Tian | Energy dissipation in twodimensional unsteady plunging breakers and an eddy viscosity model[END_REF][START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF][START_REF] Erinin | Spray generation by a plunging breaker[END_REF]. Studies such as these have identified the main controlling parameters of breaking waves, namely the breaking speed and the wave slope at breaking S = ak, where a is the wave amplitude and k the wavenumber. The bandwidth of the wave packet is also important, and the detailed kinematics before breaking, in particular a significant slowdown of the wave crest, have been discussed in order to propose breaking threshold criteria [START_REF] Banner | Linking reduced breaking crest speeds to unsteady nonlinear water wave group behavior[END_REF]Saket et al. 2017b;Pizzo & Melville 2019;[START_REF] Derakhti | A unified breaking onset criterion for surface gravity water waves in arbitrary depth[END_REF][START_REF] Fedele | Crest speeds of unsteady surface water waves[END_REF]), although we will neglect its influence from hereon in.

It follows that DNS of breaking waves can be framed in terms of a set of non-dimensional numbers. The relevant parameters are the air-water density and viscosity ratio, the wave speed and wavenumber, and amplitude. These define a wave Reynolds number, and the wave slope where λ 0 = 2π/k is the wavelength and ν is the kinematic viscosity of the water.

Similarly to turbulent DNS, numerical simulations of breaking waves are typically confined to the highest Re accessible to available computation effort, which has grown over time. [START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF]; [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]; De [START_REF] Vita | Breaking of modulated wave groups: kinematics and energy dissipation processes[END_REF] have typically used Re = 40 × 10 3 .

To consider bubble and droplet generation, the Bond number is needed,

Bo = ∆ρg σk 2 . (1.2)
∆ρ is the density difference between air and water and σ is the surface tension. Bo corresponds to the ratio between the wavelength and the capillary length scale. [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] used the Bond number to compare the numerical wavelength to experimental results. [START_REF] Deike | Capillary effects on wave breaking[END_REF] describes the wave patterns for a large range of Bo and S, discussing the energetics of parasitic capillary waves, spilling breakers and plunging breakers. As discussed in [START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF]; [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF], the breaking waves in laboratory would approach Re = 10 6 . Despite this difference in Re, DNS [START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF][START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] and LES [START_REF] Derakhti | Bubble entrainment and liquid-bubble interaction under unsteady breaking waves[END_REF][START_REF] Mcwilliams | Submesoscale currents in the ocean[END_REF] found good agreement between experiments and simulations for the non-dimensional energy dissipation due to breaking as a function of the breaker slope (see figure 1)). Nevertheless, an outstanding challenge in direct numerical simulation (DNS) is the correct numerical resolution of processes whose separation of scales increases with Re and Bo. For such simulations to capture the physics of breaking waves correctly, they must resolve all scales between and including those of energy dissipation and the formation and breakup of bubbles and droplets in a two-phase turbulent environment. This requires capturing the full physics of the problem, while retaining a qualitatively faithful representation of the breaking process in comparison with experiment. These very difficult challenges have historically limited the scope of DNS investigations, the details of whose approaches are discussed in more detail below.

Both the wave Reynolds number and the Bond number characterize the overall scale of the wave through its wavelength and phase speed, compared with viscous and capillary effects. Once the wave breaks, the turbulence it generates is controlled by the breaking slope together with the speed of the breaker, and is itself characterized by a turbulence Reynolds number, typically defined using the Taylor micro-scale Re λ , with [START_REF] Drazen | Turbulence and mixing in unsteady breaking surface waves[END_REF] typically finding values around Re λ ≈ 500. Similarly, the fragmentation processes and generation of drops and bubbles in a turbulent flow is usually analysed in terms of a Weber number, comparing the inertial stresses due to the turbulence to the surface tension.

Energetics and dimensionality of breaking waves

Breaking waves dissipate energy, generating a turbulent two-phase flow, with properties that can be related to the local breaking properties [START_REF] Duncan | An experimental investigation of breaking waves produced by a towed hydrofoil[END_REF]. The local turbulent dissipation rate due to breaking can be described by an inertial scaling [START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF],

ε = gh 3 /h (1.3)
where h is the breaking height, here consistently defined as half the distance between wave crest and trough, √ gh the ballistic velocity of the plunging breaker where g is the acceleration due to gravity. The turbulence is confined to a volume V 0 = AL c , of cross section which is generally assumed to be A πh 2 /4 [START_REF] Duncan | An experimental investigation of breaking waves produced by a towed hydrofoil[END_REF][START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF], and length of breaking crest L c , leading to an integrated dissipation rate per unit length of breaking crest, l = ρAε.

(1.4)

This scaling can be related to the initial slope, bandwidth, and speed of the wave packet in controlled laboratory experiments [START_REF] Duncan | An experimental investigation of breaking waves produced by a towed hydrofoil[END_REF][START_REF] Rapp | Laboratory measurements of deepwater breaking waves[END_REF][START_REF] Banner | Wave breaking onset and strength for two-dimensional deep-water wave groups[END_REF][START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF][START_REF] Tian | Energy dissipation in twodimensional unsteady plunging breakers and an eddy viscosity model[END_REF][START_REF] Grare | Growth and dissipation of windforced, deep-water waves[END_REF] and numerical simulations [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF][START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF][START_REF] Derakhti | Breaking-onset, energy and momentum flux in unsteady focused wave packets[END_REF]. The breaking parameter b is a non-dimensional measure of the dissipation that was introduced by [START_REF] Duncan | An experimental investigation of breaking waves produced by a towed hydrofoil[END_REF]; [START_REF] Phillips | Spectral and statistical properties of the equilibrium range in windgenerated gravity waves[END_REF], and relates to l l = bρc 5 /g, et al. 2010; [START_REF] Grare | Growth and dissipation of windforced, deep-water waves[END_REF]. Several numerical studies have confirmed this scaling and validated their approaches against this result [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF][START_REF] Deike | Lagrangian transport by breaking surface waves[END_REF][START_REF] Derakhti | Bubble entrainment and liquid-bubble interaction under unsteady breaking waves[END_REF][START_REF] Mcwilliams | Submesoscale currents in the ocean[END_REF][START_REF] Vita | Breaking of modulated wave groups: kinematics and energy dissipation processes[END_REF]. Figure 1 shows b as a function of S for a variety of experimental and numerical data, including from the present study. We note that experimental work using the linear focusing technique typically considers the linearly predicted wave slope, summed over all components, while numerical work using compact wave initialization have considered the initial slope. In all cases, the slope being used is proportional to the breaking slope, as discussed in [START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF] for experimental data and [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] for numerical data, which allows comparison between the experimental and numerical work. The differences in definitions and estimations may therefore be responsible for some of the scatter in figure 1 between the various data sets, and uncertainties in the fitting coefficients are indicated by the shaded area. Note that the scaling b ∝ S 5/2 is observed at high slopes for both the experiments and DNS. Moreover, the proportion of energy dissipated by breaking for a given slope in is similar between experiments and simulations. This fundamental model for the turbulent dissipation rate has been successfully used as the physical basis of larger-scale spectral wave models [START_REF] Romero | Spectral energy dissipation due to surface wave breaking[END_REF][START_REF] Romero | Distribution of surface wave breaking fronts[END_REF]. Moreover, we recently proposed an extension of the inertial argument to certain types of shallow-water breakers [START_REF] Mostert | Inertial energy dissipation in shallow-water breaking waves[END_REF]. It remains to determine the particular transition characteristics of the fully three-dimensional (3D) flow, and to investigate the dependence of these characteristics on the flow Reynolds number, as well as on the evolution of the ingested bubble plume. Furthermore, even aside from limitations on the maximum values of Re, Bo attainable in computation, many numerical studies have investigated two-dimensional (2D) breakers as computationally feasible analogues for the full 3D processes [START_REF] Song | A numerical study of breaking waves[END_REF][START_REF] Hendrickson | Navier-stokes simulations of unsteady small-scale breaking waves at a coupled air-water interface[END_REF][START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF][START_REF] Deike | Capillary effects on wave breaking[END_REF]. Surprisingly, despite the essentially 3D nature of the turbu-lence resulting from the breaking process, 2D breakers at the tested conditions provided a reasonable estimate of the dissipation rates for 3D breakers obtained from computation and experiment, with discrepancies sometimes as small as 5% [START_REF] Lubin | Three-dimensional large eddy simulation of air entrainment under plunging breaking waves[END_REF][START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF]. Favourable comparison with semi-empirical models as discussed above also suggests the usefulness of 2D computations for the dissipation rate [START_REF] Deike | Capillary effects on wave breaking[END_REF]. Nonetheless, the details of the 2D-3D transition physics in breaking waves constitute an open question. The present study will go some way to addressing these questions, with suggestion of a possible transition to turbulence with an associated turbulence Reynolds number.

Bubble size distributions in breaking waves

A breaking wave entrains air, which is characterized by a broad size distribution of bubbles. Direct investigation of the bubble distribution, obviously not available within a 2D study, is important to inform subgrid scale models used in LES [START_REF] Shi | Modeling quiescent phase transport of air bubbles induced by breaking waves[END_REF][START_REF] Liang | Modeling bubbles and dissolved gases in the ocean[END_REF][START_REF] Liang | Large eddy simulation of the bubbly ocean: New insights on subsurface bubble distribution and bubble-mediated gas transfer[END_REF][START_REF] Derakhti | Bubble entrainment and liquid-bubble interaction under unsteady breaking waves[END_REF] and gas transfer models [START_REF] Liang | Modeling bubbles and dissolved gases in the ocean[END_REF][START_REF] Deike | Gas transfer by breaking waves[END_REF]. [START_REF] Garrett | The connection between bubble size spectra and energy dissipation rates in the upper ocean[END_REF] proposed a turbulent break-up cascade model for the size distribution per unit volume N (r), where r is the bubble radius, as a function of the local dissipation rate ε with constant volumetric air flow-rate Q, with a dimensional analysis yielding,

N (r) ∝ Qε -1/3 r -10/3 . (1.7)
We note that a time-averaged dissipation rate ε over the breaking time has been considered when analyzing and scaling various data sets in [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF]; [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]. The corresponding breakup model assumes a turbulent inertial subrange with a direct cascade, with large bubbles injected at one end of the cascade by a notional entrainment process and turbulent fluctuations then breaking these into smaller bubbles. The lower end of the cascade is set by the Hinze scale [START_REF] Hinze | Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes[END_REF][START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF][START_REF] Perrard | Bubble deformation by a turbulent flow[END_REF],

r H = C 0 σ ρ 3/5 ε-2/5 . (1.8)
Here C 0 0.4 [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF] is a dimensionless constant. Its value is related to the critical Weber number defining bubble break-up which ranges typically from 1 to 5 [START_REF] Risso | Oscillations and breakup of a bubble immersed in a turbulent field[END_REF][START_REF] Martinez-Bazan | On the breakup of an air bubble injected into a fully developed turbulent flow. part 1. breakup frequency[END_REF][START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF][START_REF] Vejražka | Experiments on breakup of bubbles in a turbulent flow[END_REF][START_REF] Perrard | Bubble deformation by a turbulent flow[END_REF][START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF], with estimations of C 0 varying by about a factor of 2. These differences are related to variations in the experimental protocols and the large scale structure of the turbulent flow. Note also that the breaking wave problem is transient in nature, so that the Hinze scale might present variations in time, and estimations of the Hinze scale based on the averaged turbulence dissipation rate presents an added uncertainty. For all these reasons, it should be considered a soft limit. The size distribution below the Hinze scale is not addressed by [START_REF] Garrett | The connection between bubble size spectra and energy dissipation rates in the upper ocean[END_REF].

Laboratory experiments have reported measurements of the bubble size distribution under a breaking wave using various optical and acoustic techniques [START_REF] Loewen | Bubbles entrained by mechanically generated breaking waves[END_REF][START_REF] Terrill | Bubble entrainment by breaking waves and their influence on optical scattering in the upper ocean[END_REF] [START_REF] Garrett | The connection between bubble size spectra and energy dissipation rates in the upper ocean[END_REF]. Theoretical and numerical investigation has further strengthened understanding of the turbulent bubble cascade above the Hinze scale (Chan et al. 2020a,b). [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] demonstrated the ability of numerical methods to reproduce the size distribution observed experimentally and described theoretically, with an extension of the theory to constrain the mean air flow rate for increasing wave slopes. That study also noted a correspondence between the development of the entrained bubble population and the wave's energy dissipation rate. For bubbles below the Hinze scale, however, there is significant scatter between existing datasets, although [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF] suggests a relationship of ∝ r -3/2 .

The numerical studies from [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]; [START_REF] Wang | High-fidelity simulations of bubble, droplet and spray formation in breaking waves[END_REF] had limited resolution of sub-Hinze scale bubbles and were performed at Re = 40 × 10 3 , Bo = 200 with the assumption that the bubble size distributions were independent of Re, Bo, like the dissipation rate (see §1.3). The present DNS study brings to bear sophisticated methods and computational resources to test the dependence in Re, Bo, of the bubble size distribution, and to resolve the sub-Hinze bubble statistics. These constitute two of the main objectives of the present study.

Droplet size distributions in breaking waves

The mechanisms of spray generation by breaking waves have been recently reviewed by [START_REF] Veron | Ocean spray[END_REF]. Droplet size distributions have been explored experimentally in the presence of wind [START_REF] Wu | Oceanic whitecaps and sea state[END_REF][START_REF] Veron | Sea spray spume droplet production in high wind speeds[END_REF][START_REF] Ortiz-Suslow | Sea spray generation in very high winds[END_REF][START_REF] Troitskaya | The bag breakup spume droplet generation mechanism at high winds. part i: Spray generation function[END_REF] as well as for deep water breaking waves generated by linear focusing [START_REF] Erinin | Spray generation by a plunging breaker[END_REF], while numerical investigations have been made of Lagrangian transport of spume droplets in the air [START_REF] Richter | Sea surface drag and the role of spray[END_REF][START_REF] Druzhinin | The study of droplet-laden turbulent airflow over waved water surface by direct numerical simulation[END_REF][START_REF] Tang | Numerical study on the generation and transport of spume droplets in wind over breaking waves[END_REF]. However, a general theoretical model for the droplet size distribution has not been formulated.

In the context of breaking waves, spray is not created in the same manner as bubbles in the flow, being instead more analogous to atomization and fragmentation droplets [START_REF] Veron | Sea spray spume droplet production in high wind speeds[END_REF][START_REF] Troitskaya | The bag breakup spume droplet generation mechanism at high winds. part i: Spray generation function[END_REF][START_REF] Villermaux | Fragmentation versus cohesion[END_REF]. They are generated by two main mechanisms: direct ejection from wave impact and the related dynamic interface evolution, and indirect jet ejection resulting from the bursting of bubbles that were initially entrained by the breaker [START_REF] Lhuissier | Bursting bubble aerosols[END_REF]Deike et al. 2018;[START_REF] Berny | Role of all jet drops in mass transfer from bursting bubbles[END_REF]. The latter population is typically much smaller than the former [START_REF] Veron | Ocean spray[END_REF], and hence even more challenging to resolve numerically within the breaking wave event, but can be studied separately (Deike et al. 2018;[START_REF] Berny | Role of all jet drops in mass transfer from bursting bubbles[END_REF]). Separately, a major complicating factor is that spray droplet populations are typically significantly smaller than bubble populations for a given breaking wave, leading to challenges in statistical convergence of the data. For these reasons, experimental and numerical studies of droplet production by breaking waves are limited [START_REF] Wang | High-fidelity simulations of bubble, droplet and spray formation in breaking waves[END_REF][START_REF] Erinin | Spray generation by a plunging breaker[END_REF]. In this study, droplet populations are resolved over a sufficient range of length scales to allow comparison with experiment, showing good agreement in the shape of the resolved size distribution. Velocity and joint velocity-size distributions are also shown, which will aid future studies.

Outline

In this paper, we present high-resolution DNS, which mobilizes sophisticated tools and computational resources to advance the following challenges; we will show statistics spanning multiple scales of fluid behaviour for full 3D simulations which capture breaking physics as seen in laboratory experiments regarding energy dissipation, bubble and droplets size distribution. The setup is similar to [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] and is analogous to deep-water breaking waves in the laboratory obtained by focusing packets [START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF][START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF] as demonstrated by [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF], but increased resolution of the interfacial processes allows access to higher Reynolds and Bond number to describe the transition to 3D turbulence and the formation of droplets and bubbles down to scales comparable to state-of-the-art laboratory experiments. These simulations represent the current state of the art in multiphase simulations of breaking waves and further confirm that the physics of breaking waves can be profitably investigated through these high fidelity numerical data. We analyse the role of these parameters in interfacial processes, including air entrainment, bubble statistics and droplet statistics. We discuss how energy dissipation, bubble and droplet statistics seem independent of the Reynolds number above a certain value, for the strong plunging breakers, confirming the results obtained previously at lower Reynolds numbers by comparison with experimental data. Next, we investigate the role of the capillary length and other flow scales on the air entrainment and spray production, which are most likely to mediate the development of transverse instabilities in the breaking process. We emphasize that such a study is only possible thanks to improvement in adaptive mesh refinement (AMR) techniques, along with increasing computational power, which has enabled sufficiently high resolution.

The paper proceeds as follows. In §2, we describe the numerical methods and the formulation of the physical problem, the transition from the initial planar configuration to fully-developed 3D flows, and the general processes which produce entrained bubbles and ejected spray. In §3, we investigate the development of the 3D flow in direct comparisons with 2D computations as well as the role of transverse instabilities and their influence on the dissipation rate. We study the transition time and length scale of the breaking flow, from its initial 2D configuration, to final 3D turbulent one. Then, in §4, we present bubble size distribution at higher Re, Bo and numerical resolutions than those found in the numerical literature, and extend below the Hinze scale at lower Re, Bo. Droplet size and velocity distributions are presented in §5, before concluding in §6.

Problem formulation and numerical method

Basilisk library

We use the Basilisk library to solve the two-phase incompressible Navier-Stokes equations with surface tension, in 2D and 3D. The successor of the Gerris flow solver [START_REF] Popinet | Gerris: a tree-based adaptive solver for the incompressible euler equations in complex geometries[END_REF][START_REF] Popinet | An accurate adaptive solver for surface-tension-driven interfacial flows[END_REF], Basilisk is able to solve a diversity of partial differential equation systems in an adaptive mesh refinement (AMR) framework which significantly decreases the cost of high-resolution computations, allowing an efficient representation of multiscale processes. Flow advection is approximated using the Bell-Collela-Glaz method [START_REF] Bell | A second-order projection method for the incompressible navier-stokes equations[END_REF]) and the viscous terms are solved implicitly. The interface between distinct gas and liquid is described by a geometric volume-of-fluid (VOF) advection scheme, with a well-balanced surface tension treatment which mitigates the generation of parasitic currents [START_REF] Popinet | Numerical models of surface tension[END_REF]. A momentum-conserving implementation allows to avoid artefacts due to momentum "leaking" between the dense and light phases [START_REF] Fuster | An all-mach method for the simulation of bubble dynamics problems in the presence of surface tension[END_REF][START_REF] Zhang | Modeling and detailed numerical simulation of the primary breakup of a gasoline surrogate jet under non-evaporative operating conditions[END_REF]. The governing equations can be written,

∂ρ ∂t + ∇ • (ρu) = 0, (2.1) ρ ∂u ∂t + u • ∇u = -∇p + ∇ • (2µD) + ρg + σκδ s n, (2.2) ∇ • u = 0, (2.3)
where ρ, u, µ, σ, D, g are the fluid density, velocity vector, dynamic viscosity, surface tension, deformation tensor, and gravitational acceleration vector respectively. The density and viscosity are allowed to vary according to a volume fraction field c(x, t) which in these simulations takes the value of zero in the gas phase and unity in the liquid phase. The variable δ s is a Dirac delta which concentrates surface tension effects into the liquid-gas interface; κ is the curvature of the interface, and n is its unit normal vector.

Wave initialization

We consider breaking waves in deep water. The relevant physical parameters are: the liquid and gas ρ w , ρ a respectively, the respective dynamic viscosities µ w , µ a , the surface tension σ, the wavelength λ 0 , initial wave amplitude a, and gravitational acceleration g. The water depth h 0 , while finite, is assumed sufficiently large so that it does not significantly affect the breaking physics. The eight significant parameters, which are expressed in three physical dimensions, can thus be reduced into five dimensionless groups according to Buckingham's theorem; these are the density ratio ρ a /ρ w , viscosity ratio µ a /µ w , wave slope S = ak where k = 2π/λ 0 is the wavenumber, and the Bond and Reynolds numbers, as defined previously, Bo = ∆ρg/σk 2 , Re = gλ 3 0 /ν w , where ∆ρ = ρ wρ a ρ w , and ν w = µ w /ρ w is the kinematic viscosity. The wave period is T = λ 0 /c = 2π/ √ gk, where c = g/k is the linear phase speed for deep water gravity waves. The governing equations (2.1)-( 2.3) can be non-dimensionalized in terms of these groups. These definitions follow the literature, see [START_REF] Chen | Two-dimensional navier-stokes simulation of breaking waves[END_REF]; [START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF]; [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF].

The numerical resolution is indicated by the smallest cell size attained in the simulation, given by ∆ = λ 0 /2 L , where L is the maximum level of refinement used in the AMR scheme. The refinement criterion is based on both the velocity field and the VOF tracer field. The maximum resolution used in this study is L = 11, corresponding to a conventional grid of (2 11 ) 3 , or approximately 8.6 billion total cells. Under the AMR scheme, the grid size reduces to the order of 150 million cells at L = 11.

We initialize the breaking wave following [START_REF] Chen | Two-dimensional navier-stokes simulation of breaking waves[END_REF]; Iafrati (2011); [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]; [START_REF] Wang | High-fidelity simulations of bubble, droplet and spray formation in breaking waves[END_REF]; Chan et al. (2020a,b), based on an unstable third-order Stokes wave for the water velocity and zero velocity in the air. The flow is regularized in the first time step. We note that the Stokes wave solution has been derived for an irrotational, inviscid, free surface wave, hence remains an imperfect initial condition for the full two-phase flow problem, accounting for viscosity and surface tension. However, numerous studies have demonstrated that it provides an efficient and compact initialization to study the post-breaking processes. Both 2D and 3D simulations are conducted in order to investigate the transition from the laminar, planar, and essentially two-dimensional initial evolution to the final, turbulent, three-dimensional flow. Besides the dimensional difference, the 2D simulations are initialized identically to the 3D simulations. In the 3D simulations, no perturbation is used to seed the transition from planar to non-planar evolutioin of the wave; this transition is brought about by numerical noise during the breaking process.

Parameter space

The density and viscosity ratios are fixed to the values for water and air, ρ w /ρ a = 850, µ w /µ a = 51.15, and the input slope is fixed at a nominal value of S = 0.55, leaving the remaining two groups, Re, Bo to be varied. Thus, we investigate the independent effects of variation in surface tension through Bo and viscosity through Re. The fixed value of S is chosen to be sufficiently large to force the wave into a plunging breaker [START_REF] Deike | Capillary effects on wave breaking[END_REF]. We refer the reader to [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]) for an extensive study on the role of the wave slope S at constant Re,Bo. The parameters are shown in Table 1, and correspond to low (Bo=200), medium (Bo=500) and high (Bo=1000) Bond number, and low (Re=40000) and high (Re=100000) Reynolds numbers. Cases run to test grid-convergence span moderate (L = 10) and fine (L = 11) resolutions respectively. Some additional cases at a variety of Reynolds numbers are also run for the energetics comparison in §3. We reach a maximum separation of defined scales (wavelength to Hinze scale) of a factor ∼ 550. The grid size for the L = 11 case reaches 181 million cells, for a maximum runtime (excluding scheduling and queueing times) of 1.4 months and a cost of half a million CPU-hours. These highest resolution cases were run on the Stampede2 cluster at the Texas Advanced Computing Center of the University of Texas, typically on between 192 and 768 cores of the Skylake node system. (Portions of these simulations were also run on the high performance computing resources of the National Computing Center for Higher Education (CINES)). Lower resolution cases (L = 10) were run on the TigerCPU cluster at Princeton University using typically between 160 and 320 cores. Note that while these simulations are expensive, they still save several orders of magnitude over a uniform-or fixed-grid approach, which would require a prohibitively large grid size of 8.6 billion cells in the highest-resolution case.

General flow characteristics

The wave evolves in a manner similar to that seen in previous studies with similar initialization [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]. Figure 2 shows a sequence of stills at different stages of the breaking process. The initially planar wave steepens nonlinearly to a point where it locally develops a vertical interface (a,b). The wave then overturns, forming a jet which projects forward into the upstream water surface (c), and impacts onto it (d), breaking the initially planar symmetry. At this moment, a large tube of air is ingested into the liquid bulk which we refer to as the main cavity. The wave now also forms a fine-scale 3D structure at the point of impact, while ingesting the tubular cavity. This cavity persists for some time until it breaks along its length into an array of large bubbles (at t/T = 1 -1.2 e,f). In the meantime, the continuing breaking process on the surface creates a splash-up jet, as the wave proceeds into the strongly dissipative phase of the active breaking process (f) and develops into a fully-developed 3D flow (g,h) from t/T = 1.4 onwards. At late times, most of the wave energy has been dissipated in the breaking process, but the turbulent regions persist for some time, during which a very large array of spray and especially bubbles is formed (f,g,h). All the presented cases produce a large quantity of bubbles of various sizes, but spray is abundantly produced, particularly at higher Bond numbers.

These qualitative aspects of the breaking wave dynamics are crucial for a faithful representation of the breaking process. In this respect, the evolution and dynamics of the breaker closely resembles those of laboratory experiments, notwithstanding certain Bond-and Reynolds-number influences and despite the different initializations across studies. The overturning phenomenon is very similar to that seen in [START_REF] Bonmarin | Geometric properties of deep-water breaking waves[END_REF]; [START_REF] Rapp | Laboratory measurements of deepwater breaking waves[END_REF]; [START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF]; the size and shape of the main ingested cavity matches very closely that seen in a large array of theoretical, numerical and experimental studies [START_REF] Longuet-Higgins | Parametric solutions for breaking waves[END_REF][START_REF] New | A class of elliptical free-surface flows[END_REF][START_REF] New | Computations of overturning waves[END_REF][START_REF] Dommermuth | Deep-water plunging breakers: a comparison between potential theory and experiments[END_REF][START_REF] Bonmarin | Geometric properties of deep-water breaking waves[END_REF]; and the subsequent droplet-producing splash sequence closely mirrors that seen in [START_REF] Erinin | Spray generation by a plunging breaker[END_REF] (see §5). This accurate reproduction of the breaker will be further reflected in various quantitative statistical comparisons with theory and experiment in the remainder of this paper, and moreover builds high confidence in the validity of our new results.

Energetics and transition to three-dimensional turbulent flow

We determine the effect of Re (and Bo) on the development of the 3D turbulent flow underneath the breaking wave by direct comparisons of the 3D simulations with 2D counterparts.

Energy dissipation by breaking

The wave mechanical energy is E = E P + E K , where E P = V ρg(zz 0 )dV is the gravitational potential energy with a gauge z 0 chosen such that E P = 0 for the undisturbed water surface, E K = V ρ(u • u/2)dV is the kinetic energy, and the integrals are taken over the liquid volume V [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]. The instantaneous dissipation rate in the water is ε ≡ i,j ε ij where

ε ij = (ν w /2V 0 ) V (∂ i u j + ∂ j u i ) 2 dV , with ∂ i ≡ ∂/∂x i .
We decompose ε into in-and out-of-plane components ε in + ε out where ε in = i,j=x,z ε ij contains just those contributions of the deformation tensor which lie entirely in the streamwise (x) and vertical (z) directions, and ε out = ε 3Dε in comprises the remainder (i.e. the sum of terms ε iy , ε yi for i = x, y, z, y being the spanwise direction). A planar flow features only the in-plane contribution ε 3D = ε in , and a 3D flow features an additional contribution ε out (while in 2D, ε 2D ≡ ε in ).

Figure 3a shows the budget of E over time for increasing Reynolds number (Re = 10 4 , 4 × 10 4 , 10 5 ) and constant Bond number (Bo = 500), with a direct comparison between the 2D and 3D cases. For each case, E remains approximately flat at the earliest times, which correspond to the pre-broken wave where the dissipation is entirely due to the viscous boundary layer at the surface, which is properly resolved here given the high resolution in the boundary layer near the interface and has been verified for low amplitude waves (see [START_REF] Deike | Capillary effects on wave breaking[END_REF]). Breaking begins as the wave steepens and overturns at t/T 0.6 and extends through t/T = 2 and afterwards, corresponding to the impact of the wave, and the active breaking part with air entrainment and generation of turbulence. Only a small amount of energy is dissipated in the air, amounting to approximately 5% or less of the total energy budget. At small Re, viscosity is strong and the 2D and 3D budgets are in close agreement throughout the breaking process. For larger Re (smaller viscosity), the 2D and 3D curves begin to strongly diverge at a time t/T 1.2, with the discrepancy becoming more pronounced at larger Re. The percentage of energy dissipated for this high slope breaker is about 70%, close to the amount of energy dissipated in high slope plunging breakers in laboratory experiments [START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF][START_REF] Rapp | Laboratory measurements of deepwater breaking waves[END_REF].

Numerical convergence of the simulations for the energy budget and instantaneous dissipation rates are fully discussed in supplementary materials. From those results, the budget and dissipation rate at Re = 4 × 10 4 are numerically converged in 3D between L = 10, 11, as well as for Re = 10 5 between L = 10, 11 in either 2D or 3D. The comparison of dissipation rates is very good for 2D between L = 11, 12 at all Re. A 3D simulation at L = 12 at the highest Re is not feasible to run currently given the computational cost. We note that the precise time evolution of the dissipation rate is sensitive to the precise shape at impact.

Numerical resolution of characteristic dissipative scale can also be discussed. Considering Batchelor's estimate for the viscous sublayer under the prebroken wave δ ∼ λ 0 / √ Re, our results indicate that at Re = 4 × 10 4 an effective resolution of 5 cells (at L = 10) in the sublayer suffices for grid convergence. By the same estimation, we attain 6.5 cells in the viscous sublayer for Re = 10 5 at L = 11, suggesting grid convergence at this increased resolution. A resolution criteria for traditional single-phase DNS in the literature [START_REF] Pope | Turbulent flows[END_REF][START_REF] Dodd | Analysis of droplet evaporation in isotropic turbulence through droplet-resolved dns[END_REF] involves the Kolmogorov length scale η = (ν 3 w /ε) 1/4 , with k max η > 1.5 is considered sufficiently resolved, where k max = π2 L /λ 0 is the maximum resolved wavenumber. For the present simulations, for Re = 4 × 10 4 , 10 5 , at L = 11 this corresponds to k max η 3.4, 1.8 respectively, which satisfies the criterion; and is similar to resolution used in direct numerical simulations of bubble deformation in turbulence [START_REF] Farsoiya | Bubble-mediated transfer of dilute gas in turbulence[END_REF]. For details, see supplementary materials.

Without a parallel (and currently not feasible) investigation of AMR convergence with respect to uniform-grid representation at these high resolution levels, and given these are individual realizations of multiphase turbulent flows, not ensembles, some caution in the interpretation of the present data is required. Nonetheless, using these different estimate of numerical convergence, the convergence characteristics are reasonable, given the complexity of the problem.

Transition to three-dimensional turbulent flow

Figure 4 shows the time evolution of the components of the dissipation rate for increasing Re. For each case, prior to breaking, the wave is planar and ε in is the only (small) contribution to ε 3D , but the evolution of ε in , ε out on and after jet impact depends on the particular Re. For Re = 10 4 , figure 4a exhibits an almost entirely planar flow, with ε out becoming significant only late in the breaking process, when ε in , ε out both grow rapidly to their respective peak values. Before this time, the total dissipation ε 3D approximately matches ε 2D for much of the time that the flow is planar, but deviates at later times.

At higher Re = 4 × 10 4 , shown in figure 4b, 3D effects arise earlier and are much more important: ε out grows gradually from the moment of impact, and at the moment of peak dissipation, ε in and ε out are comparable. At late times, they remain similar in magnitude, suggesting that the flow has become fully 3D and turbulent by t/T = 1.3 -1.4. As before, ε 3D diverges from ε 2D at the time of rapid growth of ε in , ε out , reaching a maximum value almost double that of ε 2D .

Figure 4c, showing Re = 10 5 , is similar to figure 4b, but it does not exhibit any phase of latent planar flow where ε in ε out , and the transition to a fully 3D flow is much faster after jet impact at t/T = 0.6. Note that in this case while each of ε in , ε out is similar to ε 2D , the in-and out-of-plane contributions are not analogous to 2D processes.

Finally, figure 4d shows an overlay of each of the total instantaneous dissipation rates from figure 4a-c, suggesting that the total dissipation rate evolution and maximum value is similar between the two highest Re cases. Note however that since these cases are individual realizations of turbulent flow fields, these suggestions should be quantified further by the production and analysis of turbulent ensembles, which are prohibitively expensive to produce at these Reynolds and Bond numbers in the present investigation.

The values of ε 3D are similar for the highest Re suggesting the breaking process has achieved an asymptotic behaviour in terms of dissipation rate. The dissipation rate shown in figure 3 an asymptotic value seems to be reached at high Re, coherent with experimental estimations (grey line). Inset shows the transitions dynamics with time rescaled as (ttim)/t shear with line legend as in (a). Colors in inset are those in main plot. [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]). As such, our results are compatible with the inertial argument, experimental studies for a wide range of breakers, and previous numerical studies.

We now investigate the development to 3D flow underneath the breakers. Figure 5a shows the relative increase of the out-of-plane contributions, ε out /ε 3D with time as well as the concomitant decrease of ε in /ε 3D , for increasing Reynolds number. The terminal turbulent state is reached when either curve plateaus; this state occurs earlier for larger Re, showing the rapidity of development from planar to 3D flow. This indicates that viscosity mediates the 3D instabilities involved in the transition to turbulence at low Re.

We define a heuristic development time to 3D turbulent flow t 3D which is the time from impact until the moment that ε out /ε in = ĉ, where ĉ is some representative percentage of the turbulence dissipation rate. For ĉ = 0.5, the value of t 3D is indicated in each case of fig. 5a. The choice of ĉ is empirical. It is possible that mean velocity gradients in the streamwise-vertical plane (i.e. the in-plane mean gradients) contribute differently to the dissipation rate than out-of-plane mean gradients. On the other hand, mean gradients may in general play a much smaller role than that of turbulent fluctuations. Given these considerations, and in the absence of ensemble data with which the relative contributions of mean gradients and fluctuations can be quantified, we do not currently have a basis for prescribing a physically-informed value of ĉ. Nevertheless, we can assess the sensitivity of our time transition definition by varying ĉ across a range, which is here chosen from 0.4 to 0.6. Secondly, small fluctuations in ε out /ε 3D could affect t 3D , so we filtered the data with moving averages of window sizes of 3, 5, 7, 9, and 11 points to estimate how t 3D responds to gradual smoothing of the curve. The error bars are then estimated as the range of t 3D as estimated across both of these methods and are shown in figure 5b. We found that for each data point, the range in estimates of t 3D was determined solely by the variation of ĉ, highlighting its potential importance. The resulting error bars capture the plausible variation in transition time suggested by our data, and comment that ensemble data will shed a clearer light on this issue.

Finally, we also studied dependence of t 3D on numerical resolution for the cases Bo = 200, Re = 4 × 10 4 ; Bo = 500, Re = 2 × 10 4 ; and Bo = 500, Re = 10 5 ; we found that variation of t 3D remained within the error bars. For the case Bo = 3000, Re = 2 × 10 5 , numerical convergence cannot be assessed.

The transition to 3D turbulence can be analyzed in terms of a turbulence Reynolds number. We plot the transition time in fig. 5b for the various initial conditions as a function of two representative turbulence Reynolds number: using the integral length and velocity scales given by the breaking height h and ballistic velocity √ gh [START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF], an integral Reynolds number is Re int g 1/2 h 3/2 /ν w . The Taylor length-scale characterizing the inertial range is estimated as λ a 10/Re int and fluctuations at this scale as v = λ ε/(15ν w ) [START_REF] Sreenivasan | On the scaling of the turbulence energy dissipation rate[END_REF][START_REF] Dimotakis | Turbulent mixing[END_REF] with ε a characteristic dissipation rate taken as the peak value of ε 3D /(V 0 ρ). This yields an estimate of the turbulent Reynolds number at the Taylor micro-scale Re λ = λv/ν w 43, 69, 102 for the wave Reynolds numbers 1, 4, 10 × 10 4 .

At Re = 4 × 10 4 , the Bo = 500, 1000 points are identical, while the case Bo = 200 shows a slightly lower value of t 3D , suggesting that, for Bond numbers above 500, surface tension does not play a significant role in the transition to 3D flow and hence that Re is the main controlling parameter of this process. The inset in fig. 5b shows the relative contributions as functions of the rescaled time (tt im )/t 3D , including the different Bo, showing good collapse between all cases. The 3D transition time can be rationalized in the Re-asymptotic limit in terms of a Kelvin-Helmholtz scaling. Considering a uniform density shear layer driven by the breaker speed ≈ c = g/k over the depth of the turbulent cloud ≈ h, we get t shear 1/s where s = k KH U with U Ac and A is an O(1) constant and k KH 2/h. This shear time t shear is used to normalize the axis in fig. 5b, and since the O(1) constant is not precisely known, we indicate t shear with a shaded zone between 1 and 2 on figure 5b.

The transition time t 3D seems to plateau at the highest Reynolds number we were able to test, Re λ ≈ 50 -100. Further support of the asymptotic regime in Re number is given by considering laboratory experiments of breaking waves [START_REF] Loewen | An experimental investigation of the collective oscillations of bubble plumes entrained by breaking waves[END_REF][START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF][START_REF] Drazen | Inertial scaling of dissipation in unsteady breaking waves[END_REF][START_REF] Rapp | Laboratory measurements of deepwater breaking waves[END_REF]), with λ 0 ∼ 1 -2m, leading to Re ≈ 10 6 , and wave slopes 0.4 -0.5 inducing a turbulent flow with Re λ ≈ 500 [START_REF] Drazen | Turbulence and mixing in unsteady breaking surface waves[END_REF]. From optical and acoustic records in these experiments, we estimate the transition time t exp 3D ≈ 0.35 ± 0.1s, consistent with t exp shear ≈ t exp 3D . The transition to Reindependent flow suggests a mixing transition Reynolds number Re λ [START_REF] Sreenivasan | On the scaling of the turbulence energy dissipation rate[END_REF][START_REF] Dimotakis | Turbulent mixing[END_REF] in the flow underneath the breaking wave, supported by the similarity of ε 3D curves in fig. 4b and4c and the possibly asymptotic behaviour in fig. 5b. This suggests that Re λ 50 -100 in the developed flow corresponds with a transition to Re-independent turbulent flow under a breaking wave, which would be consistent with observations of the mixing transition in grid-generated turbulence [START_REF] Sreenivasan | On the scaling of the turbulence energy dissipation rate[END_REF] and scalar transport in turbulence [START_REF] Pullin | A vortex-based model for the subgrid flux of a passive scalar[END_REF][START_REF] Dimotakis | Turbulent mixing[END_REF]).

Air entrainment and bubble statistics

Cavity shape at entrainment

In this section, we describe air entrainment and bubble statistics. We begin by discussing the shape of the cavity at impact, which controls the size of the main cavity and the associated maximum volume of air entrained [START_REF] Lamarre | Air entrainment and dissipation in breaking waves[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]. Studies using a fully non-linear potential flow formulation, i.e. inviscid conditions and neglecting surface tension effects, have been able to reproduce the shape of the breaking wave at impact to a high level of precision [START_REF] Dommermuth | Deep-water plunging breakers: a comparison between potential theory and experiments[END_REF], with discussion on the elliptical or parametric cubic shape of the cavity [START_REF] New | A class of elliptical free-surface flows[END_REF][START_REF] Longuet-Higgins | Parametric solutions for breaking waves[END_REF]. However, these methods do not resolve the post-impact process. [START_REF] Lamarre | Air entrainment and dissipation in breaking waves[END_REF]; [START_REF] Blenkinsopp | Void fraction measurements in breaking waves[END_REF]; [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] discuss that the maximum volume of air entrained is constrained by the length of breaking crest L c and A is the cross-sectional area of the initially ingested cavity in the breaking process. In particular, A controls the amount of entrained air initially available for subsequent breakup into a bubble size distribution. It has been assumed that the cross section area of entrained air scale as A ∝ πh 2 /4 [START_REF] Duncan | An experimental investigation of breaking waves produced by a towed hydrofoil[END_REF][START_REF] Lamarre | Air entrainment and dissipation in breaking waves[END_REF][START_REF] Blenkinsopp | Void fraction measurements in breaking waves[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF], implicitly arguing that the height of the wave is large compared to the width of the jet.

As already noted in previous work, when considering a two-phase solver able to resolve post impact, moderate Bo leads to a jet thicker than observed in the laboratory [START_REF] Chen | Two-dimensional navier-stokes simulation of breaking waves[END_REF][START_REF] Song | A numerical study of breaking waves[END_REF]. Such moderate Bond numbers were nevertheless considered in most previous studies when dealing with three-dimensional breaking waves. This is because larger Bo exhibits increased separation between the wave length and Hinze scales and thus incurs a prohibitive numerical expense if all scales are to be resolved [START_REF] Wang | High-fidelity simulations of bubble, droplet and spray formation in breaking waves[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]. Here, we use the high numerical efficiency gained through AMR and increased computing power, and are thus able to resolve breakers showing greater separation between length scales. Figure 6a shows again that as Bo increases, the wave jet becomes thinner and projects further forward ahead of the wave. When increasing the Bond number, the jet at impact appears thinner and more closely similar to those observed in laboratory experiments. It is important to remark that, by comparison of the orange and red curves in figure 6a, at Bo=500, the jet thickness is independent of Reynolds number, which confirms that jet thickening is due to capillary effects. The thicker jet can be interpreted by comparing the wave height with the capillary length. For breakers in the laboratory, h ∼ 10cm and l c = γ/ρg ∼ 3mm (the capillary length) so that h/l c ≈ 33; in the DNS for Bo=200, we have h/l c ≈ 7 which indicates the importance of capillary effects. By increasing to Bo=1000, we get to h/l c ≈ 16 which is closer to laboratory conditions (but still smaller than waves from large scale breakers in the field).

V/V 0 3D 2D C = (h -π/(k √ Bo)) 2 /h 2 (e)
We therefore propose a correction of entrained area A based on the width of the jet l j . First, the cavity shape is not truly circular but closely approximates an ellipse (alternatively a parametric cubic function [START_REF] Longuet-Higgins | Parametric solutions for breaking waves[END_REF]) with an aspect ratio of √ 3 with its major axis rotated at an angle of approximately π/4 to the horizontal [START_REF] New | A class of elliptical free-surface flows[END_REF][START_REF] New | Computations of overturning waves[END_REF]. The cavity area is then A = πp 2 /(4 √ 3) where p is the major axis of the ellipse. If we now assume that, due to the thickness of the jet, the major axis is given by p 3 1/4 (h-Kl c ), where K is a positive O(1) constant, which we set to π, we obtain A = π(hπ/(k √ Bo)) 2 /4. This retrieves the usual relation for the cavity volume in the limit Bo → ∞.

Figure 6b-d show wave profiles at the moment of impact with a superimposed ' √ 3-ellipse' rotated at π/4 and with the major axis given by our estimate, 3 1/4 (hπ/(k √ Bo)). For higher Bond numbers (500, 1000), the ellipse fits very well, suggesting that our proposed cavity scaling is appropriate at high Bond numbers. Note that for the lowest Bond number ( 200), it approximates the shape of only the very rear of the cavity. This suggests that the Bo = 200 case is qualitatively distinct from higher Bo cases, in that capillary effects are sufficiently strong to change the morphology of the plunging breaker. Nevertheless, the good fit observed at higher Bond numbers supports the conjecture by [START_REF] New | A class of elliptical free-surface flows[END_REF], further supported by [START_REF] Dommermuth | Deep-water plunging breakers: a comparison between potential theory and experiments[END_REF], that the evolution of the overturning wave is independent of the details of the interior flow. Furthermore, since the √ 3-ellipse has been frequently observed in the above-cited literature, our result also confirms that this evolution is independent of the details of the initial conditions.

This leads to the cavity correction for the entrained volume, defined as the ratio of the actual entrained cavity V over its asymptotic value at high Bo number V 0 ,

C = V V 0 = (h -π/(k √ Bo)) 2 h 2 . (4.1)
This new scaling is compared with numerical data in figure 6e and shows good agreement at high Bond numbers, with weaker agreement at lower Bond numbers as expected from figure 6b-d. Note that 2D and 3D simulations are considered in figure 6e and the cavity shape is identical, since the 3D transition of the flow takes place after impact, as discussed in §3. The cavity shape is well grid converged as shown in Supplementary Materials.

Number of bubbles

We now discuss the formation of bubbles and the time evolution of their number from impact. Numerical convergence is verified for the time evolution and time averaged bubble size distribution in the Supplementary Materials.

Figure 7a shows the total number of bubbles N as a function of time (t-t im )/T . More bubbles are produced with increasing Bo, showing an order of magnitude variation in peak bubbles produced, while the production is less sensitive to Reynolds number. The total number of bubbles begins increasing at the moment of impact and peaks at the end of the active breaking stage, between 0.75T and T after impact. Particularly at higher Bond numbers, there is an increase in production rate at 0.4T which persists until ∼ 0.75T . Similar observations are made when considering only the super-Hinze scale bubbles r > r H , as shown in panel (b). The number of super-Hinze scale bubbles is much smaller than the total count, between 20 at Bo=200 to 750 at Bo=1000.

The increase of bubble production rate at (tt im )/T = 0.4 correlates with -0.5 0.0 0.5 1.0 (tt im )/T 0 5000 10000 N (t) Bo=200, Re=40000 Bo=500, Re=40000 Bo=500, Re=100000 Bo=1000, Re=100000 In both a and b, more bubbles are observed for higher Bo number, corresponding to the larger cavity. The bubble count is similar for the two Re number tested at Bo=500. (c) Detailed count breakdown for two cases, in log-log scales, showing the number of bubbles larger than the Hinze scale, r > rH and the total count, as a function of time, measured from the moment of breaking, for various cases. A nearly linear increase in number of bubbles is observed. (d) Turbulent dissipation rate as a function of time, showing both the total dissipation rate and the out-of-plane contribution, for the same cases as (c). Maximum ε is obtained when the cavity has fully collapsed.

the breakup of the main cavity. Figure 8 shows a view of the surface from below from (tt im )/T = 0.06 to 0.48. In (a) and (b), the main cavity is mostly intact, with some minor shedding of bubbles appearing off a limb of the cavity in (b). Due to the turbulence around the cavity, it deforms and ruptures dramatically in (c), creating a large number of bubbles of many sizes. The remaining parts of the cavity then destabilize further in (d), and eventually break up entirely by 0.7T after impact. Note that a significant number of bubbles is produced before this time: figure 8a shows a snapshot of the breaker from below, where many small bubbles have been entrained at the leading edge of the breaker but well before the main cavity (visible to the rear of the wave) has begun to disintegrate. Some chains of larger bubbles are also visible near the main cavity and under the primary splash-up.

Returning to figure 7, (c) shows the breakdown between sub-and super-Hinze scale bubbles for two particular cases Bo= 200 and Bo= 500, both at Re= 40000. The total count is dominated by sub-Hinze scale bubbles. The number of bubbles increases rapidly and at a roughly constant rate from the moment of impact until 0.7T or 0.8T after impact, when it begins to decay. The increase in production between 0.4T and 0.7T is subtle (on the log-log scale); before then, the bubble production rate appears to follow a broadly linear trend (indicated by the dashed black line). Finally, figure 7d shows the energy dissipation rate during the breaking process for the same two cases as (c), similarly to figure 4. Note again that the energy dissipation rate increases rapidly from 0.4T to 0.6T after impact, along with the out-of-plane contribution. The turbulence dissipation rate (as well as its out of plane contribution) is maximum when the cavity is fully broken.

This discussion suggests two effects controlling the bubble production and resulting size distribution: i) the initial air entrainment and impact, which will control initial sub-Hinze scale production, and ii) the fragmentation process of the cavity, which depends on the cavity size and the turbulence being produced during impact.

We more closely discuss the relative roles of the initial sub-Hinze production and the later multiscale fragmentation processes of the main cavity, and examine the statistics of the bubble populations. For each case, the number N and sizes of bubbles are sampled at various times t and binned by equivalent bubble radius r into bins of size ∆r, resulting in a time-dependent size distribution N (r/r H , t/T ), where r H is the Hinze scale given by (1.8) and T is the wave period, and which has been normalized by bin size such that N (r/r H , t/T )dr N (r/r H , t/T )∆r = N (t/T ), where N (t/T ) is the total number of bubbles at time t and summation is done across all radius bins.

Figure 9 shows the contours resulting from plotting N (r/r H , t/T ) over time and In each case for (tt im )/T < 0 there are no bubbles because the wave has not broken. The moment of impact corresponds with the generation of an array of sub-Hinze scale bubbles along with a single large "bubble", visible as an isolated line on the plot, which is the main cavity (see §2.4). (Individual or small numbers of similarly-sized bubbles are visible as isolated lines on the plot.) This persists until (tt im )/T 0.4; it is illustrated by figure 8a. At (tt im )/T 0.4, the cavity destabilizes and breaks into an array of large bubbles (see figure 8c,d), which themselves break up and further populate the size distribution, so that at t/T = 0.6 -0.7 there is a broad array of large and small bubbles, with the distribution weighted towards the small bubbles. At late times, ((tt im )/T = 1 onwards), the number of large bubbles reduces as they break up or reach the surface and burst. The small bubbles remain mostly entrained in the liquid for the remainder of the simulation. For a sufficiently long simulation time, all the small bubbles would eventually rise to the surface and burst, however the resolution of these bursting events would require even higher resolution (on the individual bubble) [START_REF] Berny | Role of all jet drops in mass transfer from bursting bubbles[END_REF] and are not considered here. We note that the dynamics of entrainment of the small bubbles at impact will present similarities with the physics of air entrainment by falling jets, as discussed by [START_REF] Kiger | Air-entrainment mechanisms in plunging jets and breaking waves[END_REF]. Significant qualitative differences in the distributions between the different cases are apparent only with respect to Bond number; larger Bo corresponds to a smaller Hinze scale r H , so that the distributions are generally larger relative to r H . A clear indicator is the size in r/r H of the main cavity. Re does not affect the shape of the bubble size distribution (compare the two Bo = 500 cases, at Re = 4 × 10 4 , 10 5 , since the mean turbulent dissipation rate l (which informs r H ) is not sensitive to Re for sufficiently large Re -see §3. Having discussed qualitatively the bubble production and size distribution as a function of time, we now turn to quantitative evaluations of the size distribution and its scale dependence. We focus on time averaged distributions over the active break-up time as statistical convergence of the data in the time evolution remains challenging and would require ensemble averages (requiring substantive computing time). We aim to scale the number of bubbles in the system. Figure 6 shows that the cavity shape changes at small Bo due to capillary effects, resulting in a smaller cavity, and a smaller volume of air entrained. This is confirmed by the bubble count in figure 7.

The time evolution of the bubble size distribution can be described as an extension of the model proposed by [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF], for the super-Hinze bubble distribution, based on a turbulence-buoyancy balance

N (r, t) = B AL c 2π ε(t -∆τ ) W g r -10/3 r -2/3 m , (4.2)
where A is the cross-sectional area of the initially ingested cavity in the breaking process; L c is the length of breaking crest; ε(t -∆τ ) is the energy dissipation rate and ∆τ is the time between breaker impact and peak energy dissipation rate, which corresponds to the cavity collapse time; W ≈ h/τ is a dissipation-weighted vertical mean velocity of the bubble plume over the active breaking period; with τ the active breaking period, B is a dimensionless constant. Following [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] the timescale of the cavity collapse is evaluated as ∆τ ∼ r 2/3 m ε -1/3 , where r m is the cavity size, evaluated using the scaling of the cavity length scale, r m = hl j , i.e. at high Bo number it will be independent of the Bond number (that is, the cavity of large scale breakers does not depend on surface tension), while at moderate to low Bond number, surface tension effects become important. The cross section area A controls the amount of entrained air initially available for subsequent breakup into a bubble size distribution, which we estimate from the cavity shape (see 4.1.), so that AL c ≡ V ∝ r 2 m L c . This leads to the geometric scaling N (r) ∝ r 4/3 m , which indeed indicates that the number of bubbles will increase with the size of the cavity.

Introducing the Hinze scale as characteristic length scale, eq. 4.2 can be written in a non-dimensional form as,

N (r/r H , t/T ) = B 2π ε(t -∆τ ) W g r r H -10/3 V r 3 H r H r m 2/3 . (4.3)
As described in [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF], the factor ε(t -∆τ )/(W g) describes the time evolution while the number of bubbles is determined by the strength of the break-up process and the scale separation between the initial cavity size and the Hinze scale, (V/r 3 H )(r H /r m ) 2/3 ∝ (r m /r H ) 4/3 (L c /r H ). We note that the controlling parameter in bubble break-up is the Weber number, which defines the ratio between the inertial turbulent stresses and the surface tension. When analyzing the cavity collapse, a Weber number can be defined, based on the cavity radius of r m , which depends on Bo but, as we have suggested above, approaches a constant value h/2 for sufficiently large Bo. The cavity's Weber number is then, We m = C 1 ρε -2/3 h 5/3 σ , where C 1 is a constant, and ε is the energy dissipation rate. Since the dissipation rate ε scales with the wave height h, we obtain We m = C 1 ρgS 2 σk 2 = C 1 BoS 2 ; and we further note that the scale separation r m /r H is linked to the Weber number by rm r H ∝ (W e m ) 3/5 . This links the driving Weber number of the bubble statistics and breakup processes with the Bond number and slope of the wave.

Separate studies of bubbles and droplets break-up in turbulence have demonstrated that one can observe the N (r) ∝ r -10/3 scaling in other contexts from breaking waves [START_REF] Mukherjee | Droplet-turbulence interactions and quasi-equilibrium dynamics in turbulent emulsions[END_REF][START_REF] Soligo | Breakage, coalescence and size distribution of surfactant-laden droplets in turbulent flow[END_REF][START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF], suggesting a universal character of the break-up cascade, provided the injection size is much larger than the Hinze scale, r m r H . Numerical and experimental results have shown that the number of child bubbles formed by the break-up of a large super-Hinze bubble in turbulence follows a simple power-law scaling, expressed in terms of the bubble Weber number, N ∝ (r m /r H ) α , with α between 1 and 2 [START_REF] Vejražka | Experiments on breakup of bubbles in a turbulent flow[END_REF][START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF]), which appears compatible with our results; since from eq. 4.3, V r 3

H ( r H rm ) 2/3 ∼ Lcr 4/3 m r 7/3 H ∼ ( rm r H ) 4/3 ( Lc r H ).
Note also that the DNS from [START_REF] Rivière | Sub-hinze scale bubble production in turbulent bubble break-up[END_REF] observe a nearly linear increase of the number of bubbles during the fragmentation process at high Weber number, analogous to the behaviour observed for the cavity collapse.

We consider the time-averaged version of eq. 4.3, analogous to the equation proposed by [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF], to rescale the data onto a universal scaling,

N super (r/r H ) = B 2π ε W g r r H -10/3 V r 3 H r m r H -2/3 . (4.4)
for super-Hinze scale bubbles. The sub-Hinze scale follows a r -3/2 scaling. Since the super-and sub-Hinze distributions must be continuous at the Hinze scale, we obtain We rescale the experimental distribution by the estimated cavity volume, as [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF] report a bubble size distribution n(r) in units of number of bubbles per bin size, per unit volume. For the present comparison, we consider that all bubbles are initially contained in the cavity volume V 0 . [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF]. Experimental data is plotted over r/rH on the abscissa, and scaled according to eq. 4.4. All data collapse reasonably onto a single curve given the complexity of the problem. The sub-Hinze volume is about 6% of the total entrained volume.

N sub (r/r H ) = B 2π ε W g r r H -3/2 V r 3 H r m r H -2/3 . ( 4
The time-averaged bubble size distribution, for all Re and Bo cases, over the active breaking time t/T ∈ [0 : 1.2] are shown in figure 10, and compare with the laboratory experiments from [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF]. For all Bo numbers, the bubble size distribution follows the direct cascade scaling for super-Hinze bubbles, N (r/r H ) ∝ (r/r H ) -10/3 . We resolve up to one order-of-magnitude below the Hinze scale at L = 11, in the Bo=200 case. For all cases, within this range, the size distributions have developed a shape which is clearly less steep than the super-Hinze results, close to the r -3/2 scaling, but the transition between the two regimes is not as sharp as observed in the experimental data. Note that our simulations stop at the end of the active breaking period, and as such do not describe the late-time plume evolution and steepening of the bubble size distribution, which evolves due to both degassing and further breakup, as discussed by [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF]; [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF]; [START_REF] Gaylo | Effects of power-law entrainment on bubble fragmentation cascades[END_REF]. For Bo=200, where the numerical resolution is sufficient to allow for a discussion of the sub-Hinze scale bubbles, we observe a scaling compatible with the experimental data set from [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF], N (r/r H ) ∝ (r/r H ) -3/2 . The size distribution is normalized such that N (r/r H )d(r/r H ) = N , the total number of bubbles. The partitioning in volume of air entrained is about 94% within the super-Hinze range of scale and about 6% of the air within the sub-Hinze bubbles, similar to the discussion of [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF]. Figure 10 shows that the distribution in the super-Hinze regime between the Bo=1000 and the experimental [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF] data agree reasonably well in the super-Hinze region and suggests that the asymptotic regime in Bo observed for the cavity volume in figure 10 has been reached. All data in figure 10 are reasonably well collapsed onto a single curve including the experimental data of [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF], given the uncertainties in the measurements and estimations of the various terms in the scaling model.

Droplet statistics

Stages of droplet production

We now discuss droplet production. Although all breaking waves in this study produce some droplets, large numbers of droplets only appear at larger Bo. Figure 11 shows qualitatively some of the different production mechanisms observed in these cases. Some droplets are produced immediately on impact (figure 11a); from a secondary splash-up (11b); a sustained surface splashing in the developed breaker (11c); and some jet droplets, which are partially resolved in these simulations (11d). Numerical convergence of our data is discussed in details in Supplementary Materials.

Figure 12a shows the sizes of droplets produced by the secondary splash relative to the mesh size, suggesting that many of these droplets in particular have radii of approximately the smallest mesh size, hence the have to be considered with caution. Figure 12b shows a fragmenting jet produced later in the breaking process, with only the largest droplets exhibiting a radius of more than double the mesh size. The largest droplets appear during the sustained splashing phase (corresponding to figure 11c) and statistics for such droplets are numerically converged.

The total droplet production over time is shown in Figure 13. Fewer droplets are produced for all cases compared to the bubble count (fig. 7), and for Bo = 200, fewer than 100 droplets are produced over time, which prevents any statistical convergence of the distribution. The number of droplets produced increases with Bo number, and for Bo = 500 both Re show a similar time evolution in the number of drops, with about 200 drops at most. The Bo = 1000 case shows the largest droplet counts, with many droplets produced at early times after impact and up to 800 drops. For the higher Bond number cases, figure 13a shows two prominent peaks in the droplet production. The first, sharp peak occurs at approximately the same time for both Bo = 500, 1000, at (tt im )/T 0.2. Figure 11b shows qualitatively the flow around this time for the Bo = 1000 case: shortly after the initial impact, which produces a small amount of droplets, there is a secondary impact between the splash-up and the bulk of the wave; this causes a second splash-up which projects directly upwards from the surface and produces many droplets. For Bo = 200, while this same process occurs, surface tension is too strong to allow this secondary splash-up to generate droplets. This corresponds with the first peak in figure 13a, and explains why it only appears for large Bond numbers. The peak is sharp because the droplets are produced in a single well-defined process and they are quickly destroyed as they fall back to the surface. The second peak is broader Bond number dependency seen in the sizes of the droplets produced in the first peak; that is, increased Bond number produces more droplets (as in figure 13a) as well as larger ones.

Time-averaged distribution and comparison with Erinin et al 2019

We now seek to compare the present numerical data with experiment. For this purpose, we consider the droplet size distributions time-averaged over (tt im )/T ∈ [0.2, 1]. The experimental data presented in [START_REF] Erinin | Spray generation by a plunging breaker[END_REF] are reported as a droplet count per bin size, per unit length of breaking crest. In order to compare to the numerical data, we multiply [START_REF] Erinin | Spray generation by a plunging breaker[END_REF] by the wave tank width (1.15m), which yields an absolute number of drop distribution, per unit bin-size. We consider only the Part I data from Erinin et al. which corresponds to the earlier splashing stage which is best resolved in our data and we do not consider the later drop production stage which corresponds to jet drop production. We observe a reasonable agreement between our numerical data and those of [START_REF] Erinin | Spray generation by a plunging breaker[END_REF] in the range of drop size 0.08r d /l c to r d /l c , in terms of total number of ejected drops and scaling with radius. This observation is extremely encouraging, as we note that the breaker from Erinin is at a slightly smaller wave slope than our breakers, and a slope (or wave height, or falling jet speed) dependency is expected. Note that again, the effect of Reynolds number is small, since the two Bo = 500 cases at Re = 4 × 10 4 , 10 5 collapse well.

To compare the scale of drops being produced, we normalized the drop size by the capillary length l c , and therefore present N d (r d /l c ) as a function of r d /l c , shown in figure 15. We observe a remarkable general agreement in shape in overall number of drops, as well as range of drops produced, while the data of Erinin et al. extend to smaller droplets that we are not yet able to resolve. Some Bonumber dependency is evident in the numerical data which can be attributed to the enhanced surface tension effects that reduce the fragmentation process at low Bo-number. As with the bubble size distributions, there probably exists a high-Bo regime independent of surface tension, but this critical Bo number has not yet been identified. The expected dependency in slope also complicates the analysis. Understanding these effects requires both experimental and numerical data at various slopes. However, these open questions do not reduce the importance of having achieved direct numerical simulations of drop production by a splashing process, which are well-resolved numerically and agree reasonably well with the experimental data in terms of the range of drop size produced and their total number, despite significant differences in the details of the initialization between them. We do remark that the presence of wind, not accounted for in this study, would likely generate spume, which would affect the droplet size distributions, but would require the resolution of the turbulent boundary layer forcing the wave.

Droplet velocity statistics

Lastly, we consider the statistics of droplet velocities. Figure 16a shows a contour plot of the droplet velocities (normalized by the deep water phase speed c ph = g/k) over time. It shows that smaller velocities on the order of v ∼ c ph are prevalent throughout the breaking process, with larger droplet velocities ∼ 3c ph -4c ph appearing during the secondary splash (see figure 11b) and the sustained splashing later in the breaking period (see figure 11c). Comparison with figure 14d shows these larger velocities are attained at the same time that large droplets appear. Indeed, the joint distribution of velocities and droplet radii, during the time of the sustained splashing ((tt im )/T 0.6) shown in figure 16c suggests that the highest speeds are attained by the largest droplets, though there are not many such droplets. Large droplets may also be very slow. Most droplets are small (as confirmed by the marginal size distribution, shown in figure 16b and matching earlier figures) but these vary broadly in speed. Finally, the marginal velocity distribution is shown in figure 16d, showing a peak in droplets that have low speeds ∼ c ph with a drop-off at very small speeds and a skew toward high speeds. The distribution is not governed by a power law, unlike the size distributions, but appears to be best described by a Gamma distribution (solid line) or by a log-normal distribution (dotted), both of which have been observed in many fragmentation processes [START_REF] Ling | Spray formation in a quasiplanar gas-liquid mixing layer at moderate density ratios: a numerical closeup[END_REF][START_REF] Villermaux | Fragmentation versus cohesion[END_REF].

It should be noted that the velocities presented in figure 16 are those of all droplets in the gas phase, and therefore represent droplets at all points in their ballistic trajectories. The data therefore does not in general represent only ejection speeds per se. Nevertheless, it can be assumed that the largest droplet velocities observed in figures 16a,c are those of ejecting droplets, since no larger velocities are ever observed. Thus the fastest ejection speeds in the data are of the order of 3c ph -4c ph , and they mostly occur for droplets larger than approximately 0.15l c and up to 0.4l c -0.5l c . Complete statistical separation of the just-ejected droplets from the rest of the droplet population remains to be conducted in a future study.

Concluding remarks

We have presented high-resolution simulations of breaking waves using DNS of the two-phase Navier-Stokes equations with surface tension exhibiting transition in a multiphase environment from laminar to turbulent flow, for a wide range of Reynolds numbers. By varying Bond and Reynolds numbers at high numerical resolution, we discuss the energetics of the breaker as well as statistics for bubble and droplet populations. For the energy, we have analysed the transition to 3D flow in terms of the volume-integrated dissipation rate in the water phase, and showed a Reynolds-number dependency for values of the wave Reynolds number less than 10 5 , which corresponds to a mixing transition at a turbulence Reynolds number of Re λ 100, analogous to results in a variety of canonical single-phase turbulent flows. We characterize the transition time scale, which is associated with a shear mechanism, the horizontal breaker speed and the vertical breaker height. The result thus appears generic for highly-energetic breaking waves at high slope. The shear-layer instability mechanism driving the transition is local and is expected to be independent of the type of breaker (spilling or plunging). Other features of the energetics such as a large peak in dissipation rate during the active breaking phase can be explained in terms of the breakup of the main cavity entrained by the plunging breaker. This critically contextualizes prior observations in the literature that the energetics of numerical 2D breakers approximate those of 3D breakers [START_REF] Song | A numerical study of breaking waves[END_REF][START_REF] Iafrati | Numerical study of the effects of the breaking intensity on wave breaking flows[END_REF][START_REF] Iafrati | Energy dissipation mechanisms in wave breaking processes: spilling and highly aerated plunging breaking events[END_REF][START_REF] Deike | Capillary effects on wave breaking[END_REF].

Regarding the bubble statistics, we resolve across multiple scales extending from the main cavity to below the Hinze scale, particularly at low wave Bond numbers, and find reasonable agreement with experiments [START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF] across the full range of resolved bubble sizes. We describe capillary effects on the plunging jet and ingested cavity, and characterize an asymptotic Bond number. We extend the bubble size distribution model from [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] to account for variation due to capillary effects in the size of the main cavity ingested by the breaker, and in the subsequent fragmentation and breakup cascade of the cavity. Incorporated in the scaling, and as noted by [START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF], is the close connection between the bubble statistics and the energy dissipation rate in the bulk liquid. The scaling shows good collapse of the data and, again, good agreement with experiments.

We also present statistics on the droplet populations produced by the breakers. We find good agreement in the shape of the droplet size distributions with the recent experiments of [START_REF] Erinin | Spray generation by a plunging breaker[END_REF], although some slope and Bond number effects are present and remain to be precisely quantified. Statistics on the droplet velocities are discussed, and it is found that the fastest-ejecting droplets travel at up to four times the phase speed of the wave, and are also some of the largest droplets; these are produced during the most intense splashing periods of the breaker.

The bubble and droplet size distribution seem to be both independent of the Reynolds number, once above the critical Reynolds number identified in studying the 3D turbulence transition. Consistent results in simulations and experiments for the bubble and droplet size distributions, when scaled by the characteristic length scale of the problem, reinforces the discussion in the literature [START_REF] Deike | Capillary effects on wave breaking[END_REF][START_REF] Deike | Air entrainment and bubble statistics in breaking waves[END_REF] that the details of these breakers are essentially local in the sense that whatever the initial conditions of the breaker, the dissipative, bubble and droplet properties depend only parameters of the wave at the point of breaking and not on the pre-breaking history of the wave.

We note that the results discussed here are grid converged, thanks to the use of adaptive mesh refinement techniques, which allow an effective grid size of 2048 3 grid points. These results show the ability to resolve the mixing transition in the turbulent flow in multiphase DNS of 3D breaking waves, and paves the way for realistic direct simulations of turbulent two-phase flows.

Figure 1 .

 1 Figure 1. Breaking parameters b as a function of wave slope S. Red X: Present DNS data. Red circles: DNS data of Deike et al. (2016), blue circles are DNS data from Deike et al. (2015). Black and grey data: Experimental data due to Drazen et al. (2008); Banner & Peirson (2007); Grare et al. (2013). Solid line: b = 0.4(S -0.08) 5/2 , semi-empirical result of Romero et al. (2012). Shaded area indicates the uncertainties on the scaling for b.

Figure 2 .

 2 Figure 2. Snapshot renderings of the three-dimensional (3D) breaking wave water-air interface at different times, for the case Bo = 500, Re = 100 × 10 3 , at resolution L = 11. (a) t/T = 0.37, nonlinear steepening and initial overturning, (b) t/T = 0.56, jet formation (c,d) t/T = 0.67, 0.8, impact and ingestion of main cavity (e) t/T = 1.04, splash-up of main wave and rupture of main cavity (f, g, h) t/T = 1.2, 1.36, 1.52 continuation and slowdown of main breaking process.

Figure 3 .

 3 Figure 3. Energy budget for a breaking wave. (a) Energy budgets comparing 2D and 3D for Bo = 500 and Re = 1, 4, 10 ×10 4 . (b) Corresponding instantaneous resolved dissipation rate comparing 2D and 3D for ε for Bo = 500, Re = 4 × 10 4 . All simulations are at Level 11. Grid convergence studies are presented in Supplementary Materials.

Figure 4 .

 4 Figure 4. Resolved instantaneous dissipation rates for breaking waves, showing in-plane and out-of-plane contributions to the 3D dissipation rate, along with the corresponding 2D case. Bo = 500, (a) Re = 1 × 10 4 ; (b) Re = 4 × 10 4 ; (c) Re = 1 × 10 5 . The effective resolutions for each case are 1024 3 , 2048 3 , 2048 3 respectively. (d) Overlaid total instantaneous dissipation rates for each of the cases from (a),(b),(c), showing similar dissipation rate time-evolution, especially for the two highest Re values. Larger Re corresponds with a more rapid transition from a planar initial flow to fully-developed 3D flow. ε(t) is normalized by ε0, the turbulent dissipation rate predicted by the inertial scaling argument (eq. 1.3).

Figure 5 .

 5 Figure 5. (a) Transition to fully 3D flow, measured as the relative contribution to dissipation rate εout/ε3D (solid) and εin/ε3D (dash-dotted) as a function of time for various Re number. The transition time is estimated as εout/εin = ĉ = 0.5 and indicated as vertical dotted lines. Larger Re drives more rapid transition. (b) Transition time t3D/t shear as a function of Reint for all cases:an asymptotic value seems to be reached at high Re, coherent with experimental estimations (grey line). Inset shows the transitions dynamics with time rescaled as (ttim)/t shear with line legend as in (a). Colors in inset are those in main plot.

Figure 6 .

 6 Figure 6. (a) Profiles of the volume-of-fluid interface at the moment of wave impact, obtained from 2D simulations. Blue: Bo = 200, Re = 10 5 . Orange: Bo = 500, Re = 4 × 10 4 . Red: Bo = 500, Re = 10 5 . Purple: Bo = 1000, Re = 10 5 . L = 11 for each case. Note the dependence of cavity size on Bond number, but not on Reynolds number (compare red and orange curves). (b,c,d): Wave profiles at the moment of impact with superimposed fitting ellipses for the cases, (b) Bo = 1000, (c) Bo = 500, (d) Bo = 200, with Re = 10 5 in all cases. (e) Plot of cavity area over different Bond numbers. Dashed: Corrective scaling from (4.1).

Figure 7 .

 7 Figure 7. (a) Total number of bubbles as a function of time (ttim)/T . (b) Number of bubbles of size greater than the Hinze scale, r > rH as a function of time.In both a and b, more bubbles are observed for higher Bo number, corresponding to the larger cavity. The bubble count is similar for the two Re number tested at Bo=500. (c) Detailed count breakdown for two cases, in log-log scales, showing the number of bubbles larger than the Hinze scale, r > rH and the total count, as a function of time, measured from the moment of breaking, for various cases. A nearly linear increase in number of bubbles is observed. (d) Turbulent dissipation rate as a function of time, showing both the total dissipation rate and the out-of-plane contribution, for the same cases as (c). Maximum ε is obtained when the cavity has fully collapsed.

Figure 8 .

 8 Figure 8. Snapshots of the liquid-gas interface, from below, for the case Bo = 1000, Re = 10 5 , at times (ttim)/T of (a) 0.06, (b) 0.28, (c) 0.38, (d) 0.48. Immediately after breaking, showing the fully-resolved scales of bubble phenomena. Note in particular the very small bubbles visible at and in front of the leading edge of the breaker. At later stages the air cavity collapses and lead to a wide range of bubble sizes.

Figure 9 .

 9 Figure 9. Contours of bubble size distribution over time, (a) Bo = 200, Re = 40 × 10 3 ; (b) Bo = 500, Re = 40 × 10 3 , (c) Bo = 1000, Re = 100 × 10 3 , (d) Bo = 500, Re = 100 × 10 3 . L = 11 for each case. With increasing Bond number, the main cavity size increases compared to the Hinze scale. The bubble statistics is similar for the two Re numbers at Bo=500. Small sub-Hinze scale bubbles are produced at impact, while a broad bubble cascade occurs once the cavity collapses.

4. 3 .

 3 Bubble size distribution over the active breaking time: scalings

  Bubble size distribution over the active breaking time: comparison with laboratory experiments

Figure 10 .

 10 Figure10. Time-averaged bubble size distributions over the full time window, N (r/rH ), together with the experimental data from[START_REF] Deane | Scale dependence of bubble creation mechanisms in breaking waves[END_REF]. Experimental data is plotted over r/rH on the abscissa, and scaled according to eq. 4.4. All data collapse reasonably onto a single curve given the complexity of the problem. The sub-Hinze volume is about 6% of the total entrained volume.

Figure 11 .

 11 Figure 11. Snapshots of the liquid-gas interface at different magnifications and different times, showing different stages of droplet production, for the case Bo = 1000, Re = 10 5 . (a) (ttim)/T = 0.06, splashing produced by the initial impact at the front of the breaker; (b) (ttim)/T = 0.2, secondary splash-up shortly after impact producing a peak in droplet count, (c) (ttim)/T = 0.52, sustained droplet production later in the active breaking phase, (d) (ttim)/T = 1.1, jet droplet production at late times.

Figure 12 .

 12 Figure12. Snapshots of the liquid-gas interface for two droplet production stages, showing overlaid section of the numerical mesh, for the case Bo = 1000, Re = 10 5 . (a,b) (t-tim)/T = 0.2, production of fine droplets by the secondary splash-up. Many of these droplets are resolved to less than four mesh cells per droplet diameter, for which numerical convergence is difficult to achieve. (c,d) (ttim)/T = 1.2, jet and droplet production after bubble bursting. Note that the largest droplets exceed four mesh cells per diameter.

Figure 14 .

 14 Figure 14. Contours of droplet size distribution over time. (a) Bo = 200, Re = 40 × 10 3 ; (b) Bo = 500, Re = 40 × 10 3 ; (c) Bo = 500, Re = 100 × 10 3 ; (d) Bo = 1000, Re = 100 × 10 3 . L = 11 for all cases.

Figure 15 .

 15 Figure 15. Time-average of droplet size distributions over the time window t/T ∈ [0.2, 1.0], and experimental data from Erinin et al. (2019). Experimental and numerical data is scaled consistently.

Figure 16 .

 16 Figure 16. (a) Contour of droplet velocities for the case Bo = 1000, Re = 10 5 , plotted over time on the horizontal axis and velocity normalized by the wave phase speed on the vertical. (b) The droplet size distribution at the time (ttim)/T = 0.6 averaged over a time width of ∆t/T = 0.1. Dashed line: power law with exponent -2, as in Region I of Erinin et al. (2019). (c) The velocity distribution for the same time as (b). Lines: Fits for gamma (solid) and log-normal (dotted) distributions. (d) The joint size-velocity distribution for the same time as (b).

Table 1 .

 1 Computational

	Re	Bo	L	λ0/rH	∆/rH	∆/lc	rH /lc	Cost (CPU-h)
	40 × 10 3 40 × 10 3 40 × 10 3 40 × 10 3 100 × 10 3 500 200 200 500 500 100 × 10 3 500 100 × 10 3 1000 11 11 10 11 10 11 10 100 × 10 3 1000 10	284 143 501 251 501 251 767 384	0.139 0.279 0.489 0.489 0.245 0.484 0.375 0.738	0.043 0.087 0.069 0.137 0.069 0.137 0.097 0.194	0.312 0.311 0.280 0.281 0.280 0.284 0.259 0.263	1.75 × 10 5 3.22 × 10 4 2.58 × 10 5 3.83 × 10 4 5.26 × 10 5 6.42 × 10 4 5.56 × 10 5 8.76 × 10 4
	Total Cost (CPU-h)							1.76 × 10 6

matrix of parameter space for 3D breaking waves. The slope for each case is S = 0.63 modeling a strong plunging breaker. The column labels are as follows: Re -Reynolds number; Bo -Bond number; L -maximum level of grid refinement; λ0/rH , ratio of wavelength to Hinze scale; ∆/rH -ratio of smallest grid size to Hinze scale; ∆/lc -ratio of smallest grid size to the capillary length, defined as l 2 c = 1/(k 2 Bo) where k = 2π/λ0 is the wavenumber; rH /lc -ratio of Hinze scale to capillary length.

Acknowledgements

This work was supported by the National Science Foundation (Physical Oceanography) under Grant No. 1849762 to L.D. and the Cooperative Institute for Earth System modeling between Princeton and the Geophysical Fluid Dynamics Laboratory (GFDL) NOAA. Computations were partially performed using allocation TG-OCE180010 to W.M. from the Extreme Science and Engineering Discovery Environment (XSEDE), which is supported by NSF Grant No. ACI-1053575; and the HPC resources of CINES and TGCC under the allocations 2019-A0072B07760; 2020-A0092B07760 granted by GENCI, and from the Jean Zay Grand Challenge allocation from IDRIS. Computations were also performed on resources managed and supported by Princeton Research Computing, a consortium of groups led including the Princeton Institute for Computational Science and Engineering and the Office of Information Technology's High Performance Computing Center and Visualization Laboratory at Princeton University. We are grateful to the anonymous reviewers whose comments have helped improving the quality of the manuscript.

We declare no conflict of interest.

and occurs for all cases at around (tt im )/T = 0.5 to 0.6. The state around this time is shown qualitatively for Bo = 1000 in figure 11c. It occurs as the wave proceeds through its active breaking phase and is made up of many small-scale splashing events and the bursting of large bubbles that were ingested earlier in the process. Since this process is longer and not as well-defined in space or time, the peak in 13b is accordingly broader.

Figure 13b shows the energy dissipation rates for the same cases as in (a). In contrast to the close connection between the bubble statistics and energy dissipation rate, there is no clear correlation between the dissipation rate and the droplet production.

We now discuss the droplet statistics. The droplets are gathered and binned similarly to the bubbles, into distributions N d (r d /l c , t/T ), where l c is the capillary length. The droplet populations are strongly influenced by the strength of the breaker [START_REF] Erinin | Spray generation by a plunging breaker[END_REF], and by the impact of the (ballistic) jet, particularly at early times, suggesting to use the gravity-capillary length as the relevant length scale. Note also that the lack of clear dependence on Reynolds number in the drop production suggests that viscosity does not play a role in the drop formation process. Figure 14 shows the contour maps for the droplet size distributions for the cases Bo = 200, Re = 4 × 10 4 ; Bo = 500, Re = 4 × 10 4 ; Bo = 500, Re = 10 5 ; Bo = 1000, Re = 10 5 . These corroborate the picture drawn from figure 13; there are two main peaks of droplet production, which produce short-lived drops; the first peak is sharp and the second is broader. We also observe that these peaks, and especially the second peak, are the source of large droplets. There is a slight