The geometries, formation energies, and diffusion barriers of surface and subsurface intrinsic defects in graphite are calculated using spin-polarized density-functional theory and the generalized gradient approximation with a semiempirical van der Waals (vdW) correction for dispersion interactions. The calculated formation energies and diffusion barriers of subsurface interstitial (I) atoms deviate qualitatively and quantitatively from those of surface adatoms. The same trend is found also for subsurface and adatom clusters $({\text{I}}_{2},{\text{I}}_{3})$. In spite of the semiquantitative agreement on the optimized geometries, the formation energies and diffusion barriers of surface and subsurface vacancies (V), divacancies (VV), and intimate (I-V) Frenkel pairs differ significantly from the values for the analogous defects in the bulk of graphite. This suggests limited transferability of the bulk and subsurface defect models to the surface of graphite. These findings are rationalized in terms of the balance between the covalent and vdW interaction terms at the surface, subsurface, and bulk of graphite. Finally, pairing of individual defects (adatoms, I and V) is calculated to be energetically advantageous both on the surface and in the subsurface regions. This process is shown to either saturate residual dangling bonds or produce singlet spin states, thus contributing to the quenching of residual spin polarization from damaged graphite surfaces.