In this work, we benchmark the well-controlled and numerically accurate exponential thermal tensor renormalization group (XTRG) in the simulation of interacting spin models in two dimensions. Finite temperature introduces a finite thermal correlation length $\xi$, such that for system sizes $L \gg \xi$ finite-size calculations actually simulate the thermodynamic limit. In this paper, we focus on the square lattice Heisenberg antiferromagnet (SLH) and quantum Ising models (QIM) on open and cylindrical geometries up to width W=10. We explore various one-dimensional mapping paths in the matrix product operator (MPO) representation, whose performance is clearly shown to be geometry dependent. We benchmark against quantum Monte Carlo (QMC) data, yet also the series-expansion thermal tensor network results. Thermal properties including the internal energy, specific heat, and spin structure factors, etc. are computed with high precision, obtaining excellent agreement with QMC results. XTRG also allows us to reach remarkably low temperatures. For SLH, we obtain an energy per site $u_g^$≃−0.6694(4) and a spontaneous magnetization $m_S^\simeq$0.30(1) already consistent with the ground-state properties, which is obtained from extrapolated low-T thermal data on W≤8 cylinders and W≤10 open strips, respectively. We extract an exponential divergence versus T of the structure factor S(M), as well as the correlation length ξ, at the ordering wave vector M=(π,π), which represents the renormalized classical behavior and can be observed over a narrow but appreciable temperature window, by analyzing the finite-size data by XTRG simulations. For the QIM with a finite-temperature phase transition, we employ several thermal quantities, including the specific heat, Binder ratio, as well as the MPO entanglement to determine the critical temperature $T_c$.