The W-183 nuclear shielding in a variety of tungsten polyoxometalates (POM) (Lindqvist, Anderson, decatungstates, Keggin) of different shapes and charges has been modeled by DFT calculations that take into account relativistic effects, by means of the zeroorder regular approximation (ZORA), and solvent effects, by the conductor-like screening model (COSMO) continuum method. The charge/surface area ratio (q1A) is proposed as an indicator of the charge density to which the solvation energies of all POMs are correlated in a satisfactory way. Among the various theoretical levels tested (ZORA scalar or spin-orbit, frozen-core or all-electron basis set, geometry optimization in the gas phase or in the continuum solvent, etc.), the best results are obtained when both geometry optimization in solvent and spin-orbit shielding are included (mean absolute error of delta = 35 ppm). The quality of the computed chemical shifts depends systematically on the charge density as expressed by q/A; thus, POMs with low q/A ratios display the best agreement with experimental data. The performance of the method is such that computed values can aid the assignment of the W-183 NMR spectra of polyoxotungstates, as shown by the case of alpha-[PW11TiO40](5-), whose six signals are ranked computationally so as to almost reproduce the experimental ordering even though the signals are spaced by as little as 5 ppm.