This thesis considers the application of Monte Carlo and improved timestepping methods to stochastic Lagrangian models currently used by the Met Office in Atmospheric Dispersion modelling. These models are given in the form of stochastic differential equations (SDEs) and are used to predict the spread and transport of atmospheric pollutants, such as volcanic ash, in operational forecasting, emergency response situations and scientific research. The Met Office currently uses the Standard Monte Carlo (StMC) algorithm with the Symplectic Euler timestepping method. Although the StMC method can be fast enough when a low accuracy is required, it quickly becomes expensive for higher accuracies. Another difficulty is that the Symplectic Euler method requires small timesteps to be stable.To improve the methods currently used by the Met Office we consider the Multilevel Monte Carlo (MLMC) algorithm (Giles, 2008) which has shown great potential for reducing the asymptotic computational complexity when compared to StMC in many applications. In particular, we study both numerically and theoretically, by using modified equations analysis (Shardlow, 2006), (Zygalakis, 2011), (Mueller et al., 2015), the application of the complexity theorem (Giles, 2008) which describes in which cases MLMC can reduce the asymptotic cost rate. Together with this application, we develop and implement a new algorithm for the treatment of reflective boundary conditions while simulating dispersion in the atmospheric boundary layer which gives the correct variance decay rate required by the complexity theorem. In order to reduce the overall cost we consider improved timestepping methods based on a splitting approach. With the use of a Lyapunov function and the theory from Khasminskii (2011) and Milstein and Tretyakov (2005) we also show existence and uniqueness of solutions for our models and convergence of timestepping methods.Our theoretical results and numerical methods initially apply on one-dimensional models. Later, we extend our numerical methods to higher-dimensional models where we also study the effect of a background velocity field.