Modelling of volcanic ashfall has been attempted by volcanologists but very little work has been done by mathematicians. In this thesis we show that mathematical models can accurately describe the distribution of particulate materials that fall to the ground following an eruption. We also report on the development and analysis of mathematical models to calculate the ash concentration in the atmosphere during ashfall after eruptions. Some of these models have analytical solutions. The mathematical models reported on in this thesis not only describe the distribution of ashfall on the ground but are also able to take into account the effect of variation of wind direction with elevation. In order to model the complexity of the atmospheric flow, the atmosphere is divided into horizontal layers. Each layer moves steadily and parallel to the ground: the wind velocity components, particle settling speed and dispersion coefficients are assumed constant within each layer but may differ from layer to layer. This allows for elevation-dependent wind and turbulence profiles, as well as changing particle settling speeds, the last allowing the effects of the agglomeration of particles to be taken into account.