Solubility of gypsum in the quaternary systems CaSO4-HMSO4-H2SO4-H2O (HM = Cu, Zn, Ni, Mn) are predicted up to saturated concentrations of heavy metal sulfates and to a H2SO4 concentration of 2 m by a Pitzer thermodynamic model. Experimental solubility and water activity in the subbinary and subternary systems from the literature were used for model parametrization. Then the solubility phase diagrams for the quaternary systems were predicted directly with these obtained binary and ternary model parameters. In order to verify the reliability of the predicted results, a series of solubility measurements of gypsum in these quaternary systems have been carried out at 298.15 K and the measured results were compared with the predicted ones. It was shown that the Pitzer thermodynamic model can perfectly predict the solubilities of gypsum in the quaternary systems. Meanwhile, the newly obtained experimental data were compared with limited literature data in some of the quaternary systems; good agreement was found between them. Some application examples were given based on the predicted phase diagrams.