FEM has been used to model the seismic wave equation in a 2-D linear elastic and an elasto-plastic framework. The computational model has been validated by comparing the response with that from a published code. Parallel processing capability of the model has been exhibited by using a mesh with over 100 million dof on 1040 cores. For seismic wave propagation defined by stochastic physical parameters, the effectiveness of polynomial chaos expansion is investigated as a way to characterise the stochastic aspects of wave propagation. The methodology is illustrated with three numerical examples, namely (a) a linear elastic 2-D homogeneous medium, (b) a linear elastic 2-D layered homogeneous profile, (c) a non-linear elasto-plastic 2-D homogeneous medium. The results show a considerable difference in the deterministic and estimated mean response along with major scatter manifested through standard deviation, which proves why uncertainty-quantification is important in the analysis of seismic wave propagation.