Turbulent mixing resulting from hydrodynamic instabilities, such as Rayleigh-Taylor (RT), poses a significant challenge for numerical simulations used to design inertial confinement fusion (ICF) experiments. Due to its computational efficiency, Reynolds-Averaged Navier-Stokes (RANS) modeling is popular for experiment design, but current models often inaccurately predict mixing. In this dissertation, we investigate scalar transport modeling in RT mixing using the Macroscopic Forcing Method (MFM) (Mani & Park, 2021), a numerical tool for determining closure operators. We demonstrate the application of MFM to obtain moments of the spatio-temporally nonlocal eddy diffusivity and how these measurements can inform model development. Through studies in 2D and 3D across multiple Atwood numbers, we find that nonlocality is important in modeling RT mixing, and this importance increases with Atwood number. We show how this nonlocality can be incorporated into RANS models through the development of the k-L-F model.