Dislocation modelling of an earthquake fault is of great importance due to the fact that ground surface response may be predicted by the model. However, geological features of a fault cannot be measured exactly, and therefore these features and data involve uncertainties. This paper presents a Monte Carlo based random model of faults with finite element method incorporating split node technique to impose the effects of discontinuities. Length and orientation of the fault are selected as random parameters in the domain model, and hence geometrical uncertainties are encountered. Mean and standard deviation values, as well as probability density function of ground surface responses due to the dislocation are computed. Based on analytical and numerical calculation of dislocation, two approaches of Monte Carlo simulations are proposed. Various comparisons are examined to illustrate the capability of both methods for random simulation of faults.