Most dissolved iron in the ocean is bound to organic molecules with strong conditional stability constants, known as ligands that are found at concentrations ranging from 0.2 to more than 10 nmol L− 1. In this work we report the first mechanistic description of ligand dynamics in two three-dimensional models of ocean biogeochemistry and circulation. The model for ligands is based on the concept that ligands are produced both from organic matter remineralization and phytoplankton processes, and that they are lost through bacterial and photochemical degradation, as well as aggregation and to some extent in the process of phytoplankton uptake of ligand-bound iron. A comparison with a compilation of in-situ measurements shows that the model is able to reproduce some large-scale features of the observations, such as a decrease in ligand concentrations along the conveyor belt circulation in the deep ocean, lower surface and subsurface values in the Southern Ocean, or higher values in the mesopelagic than in the abyssal ocean. Modeling ligands prognostically (as opposed to assuming a uniform ligand concentration) leads to a more nutrient-like profile of iron that is more in accordance with data. It however, also leads to higher surface concentrations of dissolved iron and negative excess ligand L⁎ in some ocean regions. This is probably an indication that with more realistic and higher ligand concentrations near the surface, as opposed to the traditionally chosen low uniform concentration, iron modelers will have to re-evaluate their assumption of low scavenging rates for iron. Given their sensitivity to environmental conditions, spatio-temporal variations in ligand concentrations have the potential to impact primary production via changes in iron limitation.