We present high-resolution three-dimensional (3-D) direct numerical simulations of breaking waves solving for the two-phase Navier–Stokes equations. We investigate the role of the Reynolds number (Re, wave inertia relative to viscous effects) and Bond number (Bo, wave scale over the capillary length) on the energy, bubble and droplet statistics of strong plunging breakers. We explore the asymptotic regimes at high Re and Bo, and compare with laboratory breaking waves. Energetically, the breaking wave transitions from laminar to 3-D turbulent flow on a time scale that depends on the turbulent Re up to a limiting value
$Re_\lambda \sim 100$
, consistent with the mixing transition in other canonical turbulent flows. We characterize the role of capillary effects on the impacting jet and ingested main cavity shape and subsequent fragmentation process, and extend the buoyant-energetic scaling from Deike et al. (J. Fluid Mech., vol. 801, 2016, pp. 91–129) to account for the cavity shape and its scale separation from the Hinze scale,
$r_H$
. We confirm two regimes in the bubble size distribution,
$N(r/r_H)\propto (r/r_H)^{-10/3}$
for
$r>r_H$
, and
$\propto (r/r_H)^{-3/2}$
for
$r<r_H$
. Bubbles are resolved up to one order of magnitude below
$r_H$
, and we observe a good collapse of the numerical data compared to laboratory breaking waves (Deane & Stokes, Nature, vol. 418 (6900), 2002, pp. 839–844). We resolve droplet statistics at high Bo in good agreement with recent experiments (Erinin et al., Geophys. Res. Lett., vol. 46 (14), 2019, pp. 8244–8251), with a distribution shape close to
$N_d(r_d)\propto r_d^{-2}$
. 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.
