Organic phosphorus (P) compounds represent a major component of soil P in many soils and are key sources of P for microbes and plants. Solution NMR (nuclear magnetic resonance spectroscopy) is a powerful technique for characterizing organic P species. However, P-31 NMR spectra are often complicated by overlapping peaks, which hampers identification and quantification of the numerous P species present in soils. Overlap is often exacerbated by the presence of paramagnetic metal ions, even if they are in complexes with EDTA following NaOH/EDTA extraction. By removing paramagnetic impurities using a new precipitation protocol, we achieved a dramatic improvement in spectral resolution. Furthermore, the obtained reduction in line widths enabled the use of multidimensional NMR methods to resolve overlapping P-31 signals. Using the new protocol on samples from two boreal humus soils with different Fe contents, 2D H-1-P-31 correlation spectra allowed unambiguous identification of a large number of P species based on their P-31 and H-1 chemical shifts and their characteristic coupling patterns, which would not have been possible using previous protocols. This approach can be used to identify organic P species in samples from both terrestrial and aquatic environments increasing our understanding of organic P biogeochemistry.