In this paper, we employ the quantum imaginary time evolution (QITE) algorithm to obtain the ground state of a two-dimensional pure $\mathbb{Z}_2$ lattice gauge theory. Specifically, we perform the (classical) numerical simulation of the QITE in the $\mathbb{Z}_2$ lattice gauge theory using the tensor networks, and the results are compared with the density matrix renormalization group (DMRG) to evaluate the accuracy.
We find that the QITE can achieve a relative error of less than $0.1\%$ for a four-plaquette system and coupling values in a regime we study.
It is also observed that the algorithmic error in the QITE decreases with increasing the number of time steps, at least for our settings, at the cost of increasing the number of gates.
Furthermore, we show that utilizing Gauss's law constraints allows us to reduce both the measurement and gate costs significantly without sacrificing simulation accuracy.

