The minimally entangled typical thermal states (METTS) are an ensemble of pure states, equivalent to the Gibbs thermal state, designed with an efficient tensor network representation in mind. In this article, we use the projected entangled pair states (PEPS) as their representation on a two-dimensional (2D) lattice. Unlike matrix product states (MPS), which for 2D systems are limited by an exponential computational barrier in the lattice size, PEPS provides a more tractable approach. To substantiate the prowess of PEPS in modeling METTS (dubbed PEPS-METTS), we benchmark it against the purification method in the context of the 2D quantum Ising model at its critical temperature. Our analysis reveals that PEPS-METTS achieves accurate results with significantly lower bond dimensions. We further corroborate this finding in the 2D Fermi-Hubbard model. At a technical level, we introduce an efficient zipper method to obtain PEPS boundary MPS needed to compute expectation values and perform sampling. The imaginary time evolution is done with the neighborhood tensor update.