A computer-implemented method for determining the stability of a fault system subjected to tectonic constraints and fluid pressure based on the rupture envelopes corresponding to the transition between a slip and a no-slip condition given a friction law of the fault surfaces, the method comprising: (a) obtaining a representation that depicts one or more faults in a region of the earth as triangulated surfaces that comprises triangulated elements, wherein the representation is obtained based at least in part on geological field data including well bore data; (b) selecting variables from among parameters of the fault system, the parameters comprising a friction coefficient μ, the stress on the triangular element σ, the cohesion C o and a stress ratio R, (c) for a given set of values of ( μ, C o, R), determining an effective shear strain energy W e of a triangular element using the Mohr-Coulomb criterion M c = μσ n + C o, the effective shear strain energy W e being determined according to W e = a. τ eff 2, where a is a friction coefficient and τ eff = τ - M c if τ >Mc, 0 otherwise; and (d) summing the strain energy of each triangle in the triangulated surfaces to yield an effective shear strain energy according to W eff = a ∑ τ eff 2; (e) incrementally varying the values ( μ, C o, R) from a minimum to a maximum value and repeating steps (b) to (d) until an effective shear strain energy is computed for all points of a cube having axes defined by ( μ, C o, R); (f) using a marching cube graphics algorithm to extract one or more iso-surfaces of the effective shear strain energy based on the summing; (g) generating, based on the extracting one or more rupture envelope images, where an effective shear strain energy value of zero corresponds to a transition between a slip and a no slip condition; and (h) generating control signals configured to control field operations associated with a subsurface volume based on the created rupture envelopes. (a) obtaining a representation that depicts one or more faults in a region of the earth as triangulated surfaces that comprises triangulated elements, wherein the representation is obtained based at least in part on geological field data including well bore data; (b) selecting variables from among parameters of the fault system, the parameters comprising a friction coefficient μ, the stress on the triangular element σ, the cohesion C o and a stress ratio R, (c) for a given set of values of ( μ, C o, R), determining an effective shear strain energy W e of a triangular element using the Mohr-Coulomb criterion M c = μσ n + C o, the effective shear strain energy W e being determined according to W e = a. τ eff 2, where a is a friction coefficient and τ eff = τ - M c if τ >Mc, 0 otherwise; and (d) summing the strain energy of each triangle in the triangulated surfaces to yield an effective shear strain energy according to W eff = a ∑ τ eff 2; (e) incrementally varying the values ( μ, C o, R) from a minimum to a maximum value and repeating steps (b) to (d) until an effective shear strain energy is computed for all points of a cube having axes defined by ( μ, C o, R); (f) using a marching cube graphics algorithm to extract one or more iso-surfaces of the effective shear strain energy based on the summing; (g) generating, based on the extracting one or more rupture envelope images, where an effective shear strain energy value of zero corresponds to a transition between a slip and a no slip condition; and (h) generating control signals configured to control field operations associated with a subsurface volume based on the created rupture envelopes.