This paper presents a mixed-interpolation hybrid model for solving coupled thermo-hygro mechanical problems. The main objective is to achieve superior accuracy in stress, moisture, and heat-flux fields, while ensuring stability in long-term computations. To this end, several improvements are introduced at all stages of the finite element solution. The first innovation concerns the development of a robust variational formulation suitable for long-term simulations, following the approach adopted in our previous works on thermo-mechanics [1,2]. A regularization technique is introduced that enables the elimination of the independent rotation field from the mechanical part of the formulation. Based on the theoretical framework of the moisture model coupled with the mechanical and thermal fields, as proposed in [3], a regularized variational formulation is obtained that involves independent fields for displacement, temperature, moisture, stress, heat flux, and moisture flux. The second innovation addresses the corresponding discrete approximation, which ensures continuity of the stress vector components in the mechanical part, as well as continuity of the heat- and moisture-flux components across element boundaries in the thermal and hygral parts. This leads to the development of enhanced finite elements, previously proposed in [1] for two-dimensional thermo-viscoplasticity using enhanced constant strain triangle (CST) elements, and in [2] for three-dimensional thermo elasticity and thermo-damage using enhanced linear tetrahedral elements. In the present formulation, displacement, temperature, and moisture fields are approximated using standard linear shape functions, while stress vector components and the corresponding fluxes are interpolated using the lowest-order Raviart–Thomas vector spaces. This approach ensures field continuity and enables smooth propagation of mechanical, thermal and hygral waves across adjacent element boundaries. The final innovation concerns the time integration scheme for enhanced finite elements in time-dependent problems. A second-order accurate mid-point time integration scheme is employed, providing improved stability and robustness for long-term simulations.