Thermal tensor networks constitute an efficient and versatile representation for quantum lattice models at finite temperatures. By Trotter-Suzuki decomposition, one obtains a D+1 dimensional TTN for the D-dimensional quantum system and then employs efficient renormalizaton group (RG) contractions to obtain the thermodynamic properties with high precision. The linearized tensor renormalization group (LTRG) method, which can be used to contract TTN efficiently and calculate the thermodynamics, is briefly reviewed and then generalized to a bilayer form. We dub this bilayer algorithm as LTRG++ and explore its performance in both finite- and infinite-size systems, finding the numerical accuracy significantly improved compared to single-layer algorithm. Moreover, we show that the LTRG++ algorithm in an infinite-size system is in essence equivalent to transfer-matrix renormalization group method, while reformulated in a tensor network language. As an application of LTRG++, we simulate an extended fermionic Hubbard model numerically, where the phase separation phenomenon, ground-state phase diagram, as well as quantum criticality-enhanced magnetocaloric effects, are investigated.