Formation of oxygen vacancies in tetragonal ZrO2 was studied using ab initio embedded cluster calculations. The vacancy formation energies were calculated, and the activation energy of vacancy migration was determined for the bulk t-ZrO2. Clusters were surrounded with Zr4+ pseudopotentials located at crystal lattice sites and embedded in the Coulomb field of a set of point charges (+4/-2) modelling the Madelung potential. The set of point charges was constructed using a new procedure specially developed to improve the accuracy of the potential. The procedure is based on introducing additional compensating charges that eliminate the quadrupole moment of the unit cell and provides convergence to the true Madelung potential. The geometry of the Zr4+ pseudopotential environment was fixed, and the positions of all inner atoms were allowed to relax. In geometry optimization, calculations were performed by the Hartree–Fock method using the LANL2DZ basis sets and the corresponding effective core potentials for Zr atoms and 6-31G basis sets for oxygen (the oxygen basis set was also added at the vacancy position). At optimized geometries, the energies were calculated by the B3LYP method, and the basis sets of the central O atoms were replaced with 6-311G*. The LANL1 effective core potentials without basis sets were used for the Zr4+ pseudopotential environment. The energy of oxygen vacancy formation in bulk t-ZrO2 was found to be 8.8 eV for the process yielding a free O atom and 6.1 eV for the process yielding an O2 molecule. The barrier for oxygen vacancy migration in bulk t-ZrO2 was found to be 1.9 eV. The formation energy of oxygen vacancy at the (001) and (101) surfaces of t-ZrO2 was 8.3 eV.