Mathematical regularization of the nonlinear terms in the Navier-Stokes equations is found to provide a systematic approach to deriving subgrid closures for numerical simulations of turbulent flow. By construction, these subgrid closures imply existence and uniqueness of strong solutions to the corresponding modelled system of equations. We will consider the large-eddy interpretation of two such mathematical regularization principles, i.e. Leray and LANS-$\alpha$ regularization. The Leray principle introduces a smoothed transport velocity as part of the regularized convective nonlinearity. The LANS-$\alpha$ principle extends the Leray formulation in a natural way in which a filtered Kelvin circulation theorem, incorporating the smoothed transport velocity, is explicitly satisfied. These regularization principles give rise to implied subgrid closures which are implemented in large eddy simulations of turbulent mixing. Comparison with filtered direct numerical simulation data and with predictions obtained from popular dynamic eddy-viscosity modelling shows that these mathematical regularization models provide considerably more accuracy at a lower computational cost than the dynamic approaches. In particular, the regularization models perform especially well in capturing the flow features characteristic of the smaller resolved scales. Variations in spatial resolution and Reynolds number establish that the Leray model is more robust but also slightly less accurate than the LANS-$\alpha$ model. The LANS-$\alpha$ model retains more of the small-scale variability in the resolved solution. However, this requires a corresponding increase in the required spatial resolution. When using second-order finite volume discretization, the potential accuracy of the implied LANS-$\alpha$ model is found to be realized by using a grid spacing that is not larger than the length scale a that appears in the definition of this model.