  The statistics of breaking wave fields is characterised within a novel
multi-layer framework, which generalises the single-layer Saint-Venant system
into a multi-layer and non-hydrostatic formulation of the Navier-Stokes
equations. We simulate an ensemble of phase-resolved surface wave fields in
physical space, where strong non-linearities including wave breaking are
modelled, without surface overturning. We extract the kinematics of wave
breaking by identifying breaking fronts and their speed, for freely evolving
wave fields initialised with typical wind wave spectra. The $\Lambda(c)$
distribution, defined as the length of breaking fronts (per unit area) moving
with speed $c$ to $c+dc$ following Phillips 1985, is reported for a broad range
of conditions. We recover the $\Lambda(c) \propto c^{-6}$ scaling without any
explicit wind forcing for steep enough wave fields. A scaling of $\Lambda(c)$
based solely on the mean square slope and peak wave phase speed is shown to
describe the modelled breaking distributions well. The modelled breaking
distributions are found to be in good agreement with field measurements and the
proposed scaling is consistent with previous empirical formulations. The
present work paves the way for simulations of the turbulent upper ocean
directly coupled with realistic breaking waves dynamics, including Langmuir
turbulence, and other sub-mesoscale processes.
