³¹P solid state NMR studies combined with DFT calculations were conducted over a chosen series of rare earth element phosphates (REEPO4s), selected on the basis of the size and magnetic properties of REEs (La, Sm, Lu and Yb). PXRD analysis revealed the presence of rhabdophane (La, Sm), monazite (La) and xenotime (Lu, Yb) phases of these phosphate compounds. The direct excitation and cross-polarization ³¹P NMR studies together with calculations confirmed the PXRD results for the abovementioned bulk structures, but also revealed presence of several local phosphorus environments on surfaces. NMR is sensitive to the atomic level local interactions, and we were able to show that the combination of experimental and theoretical NMR methods can provide information unavailable with other methods. Due to the distinct coordination of the water molecules to crystal surfaces with different Miller plane cleavages, we were able to identify from the NMR spectra the surface structures of the studied minerals. This adds to the knowledge of the bulk structures of REE phosphates and provides preliminary data for studies on coordination of various ligands on REE phosphate surfaces. This combination of experimental and computational methods can further be used for studies on surface chemistry, important for applications in catalysis and extraction of REEs from the minerals.