Multiphase flows involving fluid-fluid interfaces with surface tension underlie many processes of scientific and industrial importance. A major challenge associated with the numerical modeling of these flows comes from their multi-scale behavior, where coexisting drops, bubbles, and films can vary in size by several orders of magnitude, requiring computationally prohibitive mesh resolutions. Furthermore, current methods for representing fluid interfaces rely on mesh-size-dependent numerical errors to perform topology changes, such as breakup and coalescence. This dissertation presents various numerical methods to drastically increase the accuracy of large-scale interfacial flow simulations performed with relatively coarse mesh resolutions and presents several models to account for subgrid-scale interfacial physics. Focus will be directed towards applications in spray atomization, where a large liquid structure fragments into many smaller droplets. First, the accurate calculation of surface tension forces is addressed through advances in methods for the estimation of interfacial curvature. To decouple the breakup of liquid films from the underlying mesh size, a novel two-plane method for the representation of subgrid-thickness films is then discussed. Next, a modeling framework is proposed to predict the formation of small droplets from the breakup of these liquid films, including their diameters and initial velocities. Finally, the proposed framework is validated with a canonical drop breakup problem, where the models produce drop sizes and velocities in quantitative agreement with experimental results. The proposed methods will facilitate the efficient subgrid-scale modeling of spray formation for engineering applications.