Breaking waves entrain gas beneath the surfaces of oceans. The wave-breaking process energizes turbulent fluctuations that break bubbles in quick succession to generate a wide range of bubble sizes. Understanding this generation mechanism is of practical importance, as it contributes to a better understanding of transport processes near the ocean surface, as well as interactions of the ocean surface with solar radiation and acoustic waves. In addition, it paves the way towards the development of predictive models that will reduce the computational cost of large-scale maritime and climate simulations. It has been suggested that super-Hinze-scale turbulent breakup transfers entrained gas from large to small bubble sizes in the manner of a cascade. A theoretical basis is provided for this bubble-mass cascade through an appeal to how energy is transferred from large to small scales in the energy cascade central to single-phase turbulence theories. Ensembles of numerical simulations of breaking waves are then performed using a geometric volume-of-fluid two-phase flow solver, and used to generate the bubble statistics necessary to investigate this theoretical framework and confirm the associated theoretical findings. Reliable extraction of bubble statistics from numerical simulations requires accurate and robust identification and tracking algorithms for the dispersed phase. The identification of individual bubbles and drops traditionally relies on an algorithm used to identify connected regions. This traditional algorithm can be sensitive to the presence of spurious structures. A cost-effective refinement is proposed to maximize volume accuracy while minimizing the identification of spurious bubbles and drops. An accurate identification scheme is crucial for distinguishing bubble and drop pairs with large size ratios. The identified bubbles and drops need to be tracked in time to obtain breakup and coalescence statistics that characterize the evolution of the size distribution, including breakup and coalescence frequencies, and the probability distributions of parent and child bubble and drop sizes. An algorithm based on mass conservation is proposed to construct bubble and drop lineages using simulation snapshots that are not necessarily from consecutive time steps. These lineages are then used to detect breakup and coalescence events, and obtain the desired statistics. Accurate identification of large-size-ratio bubble and drop pairs enables accurate detection of breakup and coalescence events over a large size range. Accurate detection of successive breakup and coalescence events requires that the snapshot interval be an order of magnitude smaller than the characteristic breakup and coalescence times to capture these successive events while minimizing the identification of repeated confounding events. Together, these algorithms serve as a toolbox for detailed analysis of two-phase simulations, and enable insights into the mechanisms behind bubble and drop formation and evolution in flows of practical importance. A bubble breakup cascade requires that breakup events predominantly transfer bubble mass from a certain bubble size to a slightly smaller size on average. This property is called locality. Locality is analytically quantified by extending the population balance equation in conservative form to derive the bubble-mass transfer rate from large to small sizes. Measures of locality are proposed and used to show that scalings relevant to turbulent bubbly flows, including those previously postulated and observed in breaking-wave experiments and simulations, are consistent with a strongly local transfer rate in bubble-size space, where the influence of nonlocal contributions decays in a power-law fashion at large and small bubble sizes. These theoretical predictions reveal key physical aspects of the bubble breakup cascade phenomenology, which are crucial for the generalizability of subgrid-scale models to a variety of turbulent bubbly flows. The aforementioned algorithms allow the direct measurement of locality and thus the direct verification of these theoretical predictions via numerical simulations. Using the algorithms and analytical tools described above, relevant bubble statistics are respectively measured and analyzed in the aforementioned breaking-wave simulations as ensemble-averaged functions of time. The large-scale breakup dynamics are seen to be statistically unsteady, and two intervals with distinct characteristics are identified. In the first interval, the dissipation rate and bubble-mass flux are quasi-steady (from the point of view of the small and intermediate scales), and the theoretical analysis described above is supported by all observed statistics, including the expected -10/3 power-law exponent for the super-Hinze-scale size distribution. Strong locality is observed in the corresponding bubble-mass flux in bubble-size space, supporting the presence of a super-Hinze-scale breakup cascade. In the second interval, the dissipation rate decays, and the bubble-mass flux increases as small- and intermediate-sized bubbles become more populous. This flux remains strongly local with cascade-like behavior, but the dominant power-law exponent for the size distribution increases to -8/3 as small bubbles are also depleted more quickly. This suggests the emergence of different physical mechanisms during different phases of the breaking-wave evolution, although size-local breakup remains a dominant theme. Locality implies the presence of cascade-like behavior and supports the universality of turbulent small-bubble breakup across various turbulent bubbly flows, which simplifies the development of cascade-based subgrid-scale models to predict, for example, oceanic small-bubble statistics of practical importance. In summary, this thesis presents a toolkit for population balance analysis in two-phase flows encompassing theory and simulations, as well as algorithms to bridge the two, using turbulent bubble breakup in breaking waves as its core case study, and with an eye towards model development for subgrid structures in large-eddy simulations of turbulent two-phase flows.