A complete theory for paramagnetic relaxation enhancement (PRE) and its dependence on the magnetic field is developed for systems with electron spin S = 1, 3/2, 2, 5/2, 3 and 7/2, characterised by the presence of zero-field splitting (ZFS). The electron spin interacts through dipoledipole coupling with the nuclear spin residing in the paramagnetic complex (the inner-sphere case) as well as outside of it (the outer-sphere case). The earlier theory for S = 1 and inner-sphere interaction only is included as a special case of the present, more general approach. The theory assumes a slow reorientation of the paramagnetic complex and a lack of correlation between the rotation and translation of the complex and the electron spin dynamics. The electron spin energy level structure is determined by a combination of the Zeeman interaction and the static ZFS, and depends thus on the orientation of the complex in the magnetic field. The electron spin relaxation is described by a Redfield formulation, using the pseudorotation model for the modulation of the transient zero-field splitting. Illustrative calculations are presented, showing that the typical field dependences of the inner- and outer-sphere relaxation enhancements are in general different. The static ZFS is demonstrated to influence the magnetic field dependences by affecting the frequencies occurring in the expressions for spectral densities. The model is applied to interpret the PRE of a slowly-rotating gadolinium(iii) complex, a potential magnetic resonance imaging (MRI) contrast agent.