<p>Finite element solutions of coupled thermo-hydro-mechanical (THM) systems often exhibit non-physical oscillations in pore pressure and temperature fields when using excessively small time steps, leading to a significant reduction in the accuracy, efficiency, and stability of the numerical calculations. However, the study of the oscillation mechanisms and the critical minimum time-step criteria in the coupled THM finite element analysis is still notably unexplored. This work establishes an analytical framework deriving multi-field non-oscillatory criteria for equal-order linear and quadratic interpolation schemes under both consistent and lumped mass matrices through the Discrete Maximum Principle (DMP) and Monotonicity Principle (MP). Numerical experiments are conducted to verify these theoretical results, and the characteristics of nonphysical oscillations in THM coupling problems are explored. The results show that the generation of oscillations is not only related to the time step but also to the type of load on the boundary. The field oscillations in the THM coupled system interact with each other, and the numerical oscillations are characterized by the smaller time step, the stronger and wider range of oscillations, and the numerical solutions gradually converge to a stable value along the direction away from the loads. The proposed critical minimum time-step sizes are accurate and reliable, which are affected by the element type, mesh size, mass matrix scheme, time-domain discretization scheme, and material parameters. Both theoretical and numerical analyses confirm that the linear elements have better global stability than the quadratic elements, while the lumped mass matrix has stronger oscillation suppression compared to the consistent formulations. These findings establish essential guidelines for designing robust multiphysics computational frameworks.</p>