Fault zones are structurally highly spatially heterogeneous and hence extremely complex. Observations of fluid flow through fault zones over several scales show that this structural complexity is reflected in the hydrogeological properties of faults. Information on faults at depth is scarce, hence, it is highly valuable to understand the controls on spatial and temporal fault zone development. In this paper we increase our understanding of fault damage zone development in crystalline rocks by dynamically simulating the growth of single and multiple splay fractures produced from failure on a pre-existing fault. We present a new simulation model, MOPEDZ (Modeling Of Permeability Evolution in the Damage Zone surrounding faults), that simulates fault evolution through solution of Navier's equation with a combined Mohr-Coulomb and tensile failure criteria. Simulations suggest that location, frequency, mode of failure and orientation of splay fractures are significantly affected both by the orientation of the fault with respect to the maximum principal compressive stress and the conditions of differential stress. Model predictions compare well with published field outcrop data, confirming that this model produces realistic damage zone geometries. Copyright 2007 by the American Geophysical Union.