Variational minimization of tensor network states enables the exploration of low energy states of lattice gauge theories. However, the exact numerical evaluation of high-dimensional tensor network states remains challenging in general. In [E. Zohar and J. I. Cirac, Phys. Rev. D 97, 034510 (2018)] it was shown how, by combining gauged Gaussian projected entangled pair states with a variational Monte Carlo procedure, it is possible to efficiently compute physical observables. In this paper we demonstrate how this approach can be used to investigate numerically the ground state of a lattice gauge theory. More concretely, we explicitly carry out the variational Monte Carlo procedure based on such contraction methods for a pure gauge KogutSusskind Hamiltonian with a Z3 gauge field in two spatial dimensions. This is a first proof of principle to the method, which provides an inherent way to increase the number of variational parameters and can be readily extended to systems with physical fermions.