A projected entangled pair state (PEPS) with ancillas can be evolved in imaginary time to obtain thermal states of a strongly correlated quantum system on a two-dimensional lattice. Every application of a Suzuki-Trotter gate multiplies the PEPS bond dimension D by a factor k. It has to be renormalized back to the original D. In order to preserve the accuracy of the Suzuki-Trotter (ST) decomposition, the renormalization in principle has to take into account full environment made of the new tensors with the bond dimension k D. Here, we propose a self-consistent renormalization procedure operating with the original bond dimension D, but without compromising the accuracy of the ST decomposition. The iterative procedure renormalizes the bond using full environment made of renormalized tensors with the bond dimension D. After every renormalization, the new renormalized tensors are used to update the environment, and then the renormalization is repeated again and again until convergence. As a benchmark application, we obtain thermal states of the transverse field quantum Ising model on a square lattice, both infinite and finite, evolving the system across a second-order phase transition at finite temperature.