Summary: We consider the evolution of the bulk bubble-size distribution \(N(a,t)\) of large bubbles (Weber number \(We \gg 1)\) under free-surface entrainment described generally by an entrainment size distribution \(I(a)\) with power-law slope \(\gamma\) and large-radius cutoff \(a_{\max}\). Our main focus is the interaction between turbulence-driven fragmentation and free-surface entrainment, and, for simplicity, we ignore other mechanisms such as degassing, coalescence and dissolution. Of special interest are the equilibrium bulk distribution \(N_{eq}(a)\), with local power-law slope \(\tilde{\beta }_{eq}(a)\), and the time scale \(\tau_c\) to reach this equilibrium after initiation of entrainment. For bubble radii \(a\ll a_{\max}\), we find two regimes for the dependence of \(N_{eq}(a)\) on the entrainment distribution. There is a weak injection regime for \(\gamma \geq -4\), where \(\tilde{\beta}_{eq}(a)=-10/3\) independent of the entrainment distribution; and a strong injection regime for \(\gamma <-4\), where the power-law slope depends on \(\gamma\) and is given by \(\tilde{\beta}_{eq}(a)=\gamma +2/3\). The weak regime provides a general explanation for the commonly observed \(-10/3\) power law originally proposed by C. Garrett et al. [“The connection between bubble size spectra and energy dissipation rates in the upper ocean”, J. Phys. Oceanogr. 30, No. 9, 2163–2171 (2000; doi:10.1175/1520-0485(2000)030<2163:TCBBSS>2.0.CO;2)], and suggests that different weak entrainment mechanisms can all lead to this result. For \(a\sim a_{\max}\), we find that \(N_{eq}(a)\) exhibits a steepening deviation from a power law due to fragmentation and entrainment, similar to what has been observed, but here absent other mechanisms such as degassing. The evolution of \(N(a,t)\) to \(N_{eq}(a)\) is characterised by the critical time \(\tau_c \propto C_f \varepsilon^{-1/3} a_{\max}^{2/3}\), where \(\varepsilon\) is the turbulence dissipation rate and \(C_f\) is a new constant that quantifies the dependence on the daughter size distribution in a fragmentation event. For typical breaking waves, \( \tau_c\) can be quite small, limiting the time \(t\lesssim \tau_c\) when direct measurement of \(N(a,t)\) might provide information about the underlying entrainment size distribution.