The integrated kMC-DR method is applied to simulations of zirconia film growth in the ALD process. The kMC-DR method combines conventional kinetic Monte Carlo approach with geometry optimization steps (dynamic relaxation). We found that this method is suitable for the description of the formation of non-regular structures like amorphous films. The steady state film growth rate obtained from the simulations is in a reasonable agreement with experimental results for zirconia film growth rate in an ALD reactor, thus confirming the effect of steric hindrances on the film growth rate. We also found that the zirconium coordination number increases in course of film deposition. This feature of the film deposition process was not reproduced in previous models of the ZrO₂-Si interface. The results obtained demonstrate formation of non-uniform ZrO₂ film at the first ALD cycles because of the mismatch between the oxidized silicon support and the zirconia bulk structure. Residual terminal hydroxyl groups are also formed at the ZrO₂-Si/SiO₂ interface, and we note that such groups could significantly affect band lineup at the interface in accordance with similar results for the Si-SiO₂ interface.