Dempsey et al. (this issue) have presented validation of a coupled thermal-hydrological-mechanical model through a comparative study of shear stimulation in geothermal fields at Desert Peak, Nevada, USA and Ngatamariki, New Zealand. Values of model parameters obtained from such validation studies are often used for extrapolating the model results beyond the domain of available experiments. Thus it is important to consider the sensitivity of these results to the specifics of the model setup as well as uncertainties in input parameters. In this presentation, we consider the sensitivity of history-matched parameter values to the numerical meshes on which computations are performed. During well stimulation, it has been noted that injectivity varies with time according to a power law, i.e., ࡵࡵ ∝ ࢔࢚, with ࢔ ranging between 0.3 and 0.7. Dempsey et al. propose that ࢔ depends on the geometry of the stimulated region. However, this result is dependent on permeability enhanced according to a Mohr-Coulomb failure criterion, which is informed by the stress solution. It is well known that that the finite element method, in linear elastostatics, displays optimal rates of convergence in the L2 norm of stress error with mesh refinement (Zienkiewicz and Taylor, 1994). However, the situation is significantly more complicated in fully-coupled THM modelling. Stress-induced permeability changes affect both the fluid mass balance and energy/enthalpy balance equations and the subsequent convergence of the entire coupled system of equations. Relying on solutions found on an un-converged discretization could result in significant errors. We evaluate this dependence for several grid geometries to establish the robustness of the model findings.