The present thesis deals with a number of challenges in the field of large eddy simulation (LES). These include the performance of subgrid-scale (SGS) models at fairly high Reynolds numbers and coarse resolutions, passive scalar and stochastic modeling in LES. The fully-developed turbulent channel flow is used as the test case for these investigations. The advantage of this particular test case is that highly accurate pseudo-spectral methods can be used for the discretization of the governing equations. In the absence of discretization errors, a better understanding of the subgrid-scale model performance can be achieved. Moreover, the turbulent channel flow is a challenging test case for LES, since it shares some of the common important features of all wall-bounded turbulent flows. Most commonly used eddy-viscosity-type models are suitable for moderately to highly-resolved LES cases, where the unresolved scales are approximately isotropic. However, this makes simulations of high Reynolds number wall-bounded flows computationally expensive. In contrast, explicit algebraic (EA) model takes into account the anisotropy of SGS motions and performs well in predicting the flow statistics in coarse-grid LES cases. Therefore, LES of high Reynolds number wall-bounded flows can be performed at much lower number of grid points in comparison with other models. A demonstration of the resolution requirements for the EA model in comparison with the dynamic Smagorinsky and its high-pass filtered version for a fairly high Reynolds number is given in this thesis. One of the shortcomings of the commonly used eddy diffusivity model arises from its assumption of alignment of the SGS scalar flux vector with the resolved scalar gradients. However, better SGS scalar flux models that overcome this issue are very few. Using the same methodology that led to the EA SGS stress model, a new explicit algebraic SGS scalar flux model is developed, which allows the SGS scalar fluxes to be partially independent of the resolved scalar gradient. The model predictions are verified and found to improve the scalar statistics in comparison with the eddy diffusivity model. The intermittent nature of energy transfer between the large and small scales of turbulence is often not fully taken into account in the formulation of SGS models both for velocity and scalar. Using the Langevin stochastic differential equation, the EA models are extended to incorporate random variations in their predictions which lead to a reasonable amount of backscatter of energy from the SGS to the resolved scales. The stochastic EA models improve the predictions of the SGS dissipation by decreasing its length scale and improving the shape of its probability density function.