Tensor network methods are routinely used in approximating various equilibrium and non-equilibrium scenarios, with the algorithms requiring a small bond dimension at low enough time or inverse temperature. These approaches so far lacked a rigorous mathematical justification, since existing approximations to thermal states and time evolution demand a bond dimension growing with system size. To address this problem, we construct PEPOs that approximate, for all local observables, $i)$ their thermal expectation values and $ii)$ their Heisenberg time evolution. The bond dimension required does not depend on system size, but only on the temperature or time. We also show how these can be used to approximate thermal correlation functions and expectation values in quantum quenches.