Based on bounding surface plasticity theory, both the size and shape of the bounding surface are allowed to change with temperature in the thermo-mechanical model. The thermo-bounding surface model can capture many vital aspects of thermo- mechanical behavior, including volume changes during heating and cooling, and thermal effects on shear behavior. In order to analyse thermo-fluid-solid coupling phenomena, a thermo-bounding surface model is implemented in the software COMSOL. The implementation techniques include modification of the built-in material equations of COMSOL, coupling analysis of pore fluid flow and stress, the addition of partial differential equations to compute plastic strains for any stress path within the bounding surface. The developed model was used to analyze one-dimensional consolidation, isotropic compression test at constant temperature, drained heating and cooling tests and undrained triaxial compression tests at different temperatures. By comparing with the computed and measured results, it is evident that the proposed model is capable to capture key features of the coupled soil behavior.