  A Gibbs operator $e^{-\beta H}$ for a 2D lattice system with a Hamiltonian
$H$ can be represented by a 3D tensor network, the third dimension being the
imaginary time (inverse temperature) $\beta$. Coarse-graining the network along
$\beta$ results in an accurate 2D projected entangled-pair operator (PEPO) with
a finite bond dimension. The coarse-graining is performed by a tree tensor
network of isometries that are optimized variationally to maximize the accuracy
of the PEPO. The algorithm is applied to the two-dimensional Hubbard model on
an infinite square lattice. Benchmark results are obtained that are consistent
with the best cluster dynamical mean-field theory and power series expansion in
the regime of parameters where they yield mutually consistent results.
