A computer implemented method for modeling of faulting and fracturing uses "small scale rules" to produce large-scale results. One part of the method is a user interface for inputting deformations, preexisting faults and fractures, and material rock properties. The second part of the software is the code that solves the motion of each point or node in the subsurface volume defined by the user interface. The model may be defined in one of three modes: an aerial mode, in which the model is 2-dimensional with the material and a substrate on a horizontal plane; a cross-sectional mode similar to the aerial mode except that the nodes are in a vertical cross section and gravity is included in the model; and a 3-D model that is an extension into a third dimension of the 2-D model and deformation may be applied to the bottom and four sides of the material region.A modified over-relaxation approach, wherein the over-relaxation is concentrated in those nodes where the greatest movement occurs, is used to solve for the deformation. This significantly speeds up the computation time. The model is "conditioned" to increase the likelihood that the deformation pattern resulting from the simulation of the deformation will at least duplicate an observed large-scale deformation. As an aid to the simulation, an "anticipate" step provides a quick solution to the deformation without including the effects of faulting.