Numerical modeling of advanced propulsion systems such as the Internal Combustion Engine (ICE) is of great interest to the community due to the magnitude of compute/algorithmic challenges. Fuel spray atomization, which determines the rate of fuel-air mixing, is a critical limiting process for the phenomena of combustion within ICEs. Fuel spray atomization has proven to be a formidable challenge for the state-of-the-art numerical models due to its highly transient, multi-scale, and multi-phase nature. Current models for primary atomization employ a high degree of empiricism in the form of model constants. This level of empiricism often reduces the art of predictive modeling, in the case of sprays to a mere data-fitting exercise to experimental observations by tuning model constants. In this research, first a series of full three-dimensional (3D) Computational Fluid Dynamics (CFD) studies are presented that examine the factors affecting spray atomization, including the effects of nozzle geometry, transient injector needle motion and in-cylinder thermodynamic conditions on the spray atomization behavior. Informed by this study, a novel reduced order model - Eulerian Lagrangian Mixing Oriented (ELMO), with experimentally-informed inputs is presented that aims to replace current primary atomization models such as the Kelvin-Helmholtz Rayleigh-Taylor (KHRT). The ELMO fuel spray atomization model is coupled to a full 3D gas phase solver, and validation studies on single/multi-hole injector configurations and spray conditions (diesel/gasoline) are presented. In addition, validation studies for ELMO are presented for a standard diesel engine. In the last section of this thesis, the applicability of data-driven machine learning (ML) models are explored. First, a ML based turbulence closure is presented for the four-stroke Darmstadt engine and performance against experimentally validated ground-truth data measured. Sensitivity studies indicate the data-driven model preserves the functional characteristics of the turbulence closure, consistent with theory. In the final study, a machine learning based surrogate model for estimating the discretization error of coarser meshes is developed and integrated into OpenFOAM. This surrogate model enhances the quality of coarse-mesh CFD simulations by adding necessary source terms to retrieve the 'lost' resolved scale information, thereby providing a higher fidelity solution at almost the same compute cost. This model is tested on different OpenFOAM solvers, including for an engine cold flow setting, and the generalizability of the data-driven framework commented upon.