There are two big challenges which restrict the extensive application of fully coupled geomechanics-reservoir modeling. The first challenge is computational effort. Consider a 3-D simulation combining pressure and heat diffusion, elastoplastic mechanical response, and saturation changes;​ each node has at least 5 degrees of freedom, each leading to a separate equation. Furthermore, regions of large p, T and σ′ gradients require small-scale discretization for accurate solutions, greatly increasing the number of equations. When the rock mass surrounding the reservoir region is included, it is represented by many elements or nodes. These factors mean that accurate analysis of realistic 3-D problems is challenging, and will so remain as we seek to solve larger and larger coupled problems involving nonlinear responses.

To overcome the first challenge, the displacement discontinuity method is introduced wherein a large-scale 3-D case is divided into a reservoir region where Δp, ΔT and non-linear effects are critical and analyzed using FEM, and an outside region in which the reservoir is encased where Δp and ΔT effects are inconsequential and the rock may be treated as elastic, analyzed with a 3D displacement discontinuity formulation. This scheme leads to a tremendous reduction in the degrees of freedom, yet allows for reasonably rigorous incorporation of the reactions of the surrounding rock.

The second challenge arises from some forms of numerical instability. There are actually two types of sharp gradients implied in the transient advection-diffusion problem:​ one is caused by the high Peclet numbers, the other by the sharp gradient which appears during the small time steps due to the transient solution. The way to eliminate the spurious oscillations is different when the sharp gradients are induced by the transient evolution than when they are produced by the advective terms, and existing literature focuses mainly on eliminating the spurious spatial temperature oscillations caused by advection-dominated flow.

To overcome the second challenge, numerical instability sources are addressed by introducing a new stabilized finite element method, the subgrid scale/gradient subgrid scale (SGS/GSGS) method.