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, (31)P 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 (31)P signals. Using the new protocol on samples from two boreal humus soils with different Fe contents, 2D (1)H-(31)P correlation spectra allowed unambiguous identification of a large number of P species based on their (31)P and (1)H 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.