Numerical simulations of salt cavern stress and strain under pressure fluctuation are challenging due to a combination of factors. Firstly, the solution-mined salt caverns within the heterogeneous geological layers always have quite complex geometries, motivating the use of tetrahedral mesh with local refinements. Secondly, salt rocks are known to creep under deviatoric stresses, meaning that deformations take place at constant volume (isochoric). Simulating isochoric deformations with tetrahedral meshes is particularly problematic for finite element (FE) formulations that do not satisfy the Ladyzhenskaya-Babuška-Brezzi (LBB) stability condition. This problem can indeed be avoided by employing a hexahedral mesh with the simple scheme such as B-bar method. However, hexahedral meshes are not convenient for complex cavern geometries. Alternatively, resolving the stability of the method on tetrahedral meshes is advantageous for that they are truly flexible for any cavern shape. This work presents a stabilized mixed FE formulation for low-order tetrahedral meshes that incorporates all relevant deformation mechanisms for salt rocks. The stabilization consists of enriching the displacement FE approximation in the mean stress formulation. A key step in this enrichment is to obtain an approximation for the Laplacian of the displacement field that accounts for inelastic strains. This is successfully achieved by using the Physical Influence Scheme (PIS) together with the concept of secant Young’s modulus [1], which promotes local stabilization where necessary. When combined with an accurate calculation of the scaling factor (parameter h) for the Laplacian approximation [2], this stabilization technique is shown to produce oscillation-free and physically consistent results without any tuning parameters. The proposed technique is analyzed in many relevant test cases for salt cavern simulations, and the results demonstrate its effectiveness in eliminating spurious numerical oscillations with low-order tetrahedral meshes.