  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
thermal correlation length, which justifies the analysis of finite system size
for the sake of numerical efficiency. 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 at low temperature an
energy per site $u_g^*\simeq -0.6694(4)$ and a spontaneous magnetization
$m_S^*\simeq0.30(1)$, which is already consistent with the ground state
properties. We extract an exponential divergence vs. $T$ of the structure
factor $S(M)$, as well as the correlation length $\xi$, at the ordering wave
vector $M=(\pi,\pi)$, which represents the renormalized classical behavior and
can be observed over a narrow but appreciable temperature window, by analysing
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$.
