The suitability of tensor network ansatzes for the description of physically relevant states in one dimensional lattice gauge theories (LGT) has been demonstrated in the last years by a large amount of systematic studies, including abelian and non-abelian LGTs, and including scenarios where traditional Monte Carlo approaches fail due to a sign problem. While this establishes a solid motivation to extend the program to higher dimensions, a similar systematic study in two dimensions using PEPS requires dealing with specific considerations. Besides a larger computational costs associated to the higher spatial dimension, the presence of plaquette terms in LGTs hinders the efficiency of the most up-to-date PEPS algorithms. With a newly developed update strategy, nevertheless, such terms can be treated by the most efficient techniques. We have used this method to perform the first ab initio iPEPS study of a LGT in 2+1 dimensions: a Z3 invariant model, for which we have determined the phase diagram.