The magnetic properties of paramagnetic species with spin S > 1/2 are parameterized by the familiar g tensor as well as "zero-field splitting" (ZFS) terms that break the degeneracy between spin states even in the absence of a magnetic field. In this work, we determine the mean values and distributions of the ZFS parameters D and E for six Gd(iii) complexes (S = 7/2) and critically discuss the accuracy of such determination. EPR spectra of the Gd(iii) complexes were recorded in glassy frozen solutions at 10 K or below at Q-band (∼34 GHz), W-band (∼94 GHz) and G-band (240 GHz) frequencies, and simulated with two widely used models for the form of the distributions of the ZFS parameters D and E. We find that the form of the distribution of the ZFS parameter D is bimodal, consisting roughly of two Gaussians centered at D and -D with unequal amplitudes. The extracted values of D (σD) for the six complexes are, in MHz: Gd-NO3Pic, 485 ± 20 (155 ± 37); Gd-DOTA/Gd-maleimide-DOTA, -714 ± 43 (328 ± 99); iodo-(Gd-PyMTA)/MOMethynyl-(Gd-PyMTA), 1213 ± 60 (418 ± 141); Gd-TAHA, 1361 ± 69 (457 ± 178); iodo-Gd-PCTA-[12], 1861 ± 135 (467 ± 292); and Gd-PyDTTA, 1830 ± 105 (390 ± 242). The sign of D was adjusted based on the Gaussian component with larger amplitude. We relate the extracted P(D) distributions to the structure of the individual Gd(iii) complexes by fitting them to a model that superposes the contribution to the D tensor from each coordinating atom of the ligand. Using this model, we predict D, σD, and E values for several additional Gd(iii) complexes that were not measured in this work. The results of this paper may be useful as benchmarks for the verification of quantum chemical calculations of ZFS parameters, and point the way to designing Gd(iii) complexes for particular applications and estimating their magnetic properties a priori.