Complexation with organic ligands controls the solubility of iron in the surface ocean; together with particle adsorption it determines bioavailability and residence time of dust-deposited iron. Assumptions on ligand abundance are therefore critical for biogeochemical models that include iron. Instead of fixing organic ligands to an observed mean concentration, a mechanistic description of ligand dynamics has been implemented in a 3D model. Two classes of ligands are considered: strong and weak, with different biogenic sources, photochemical reactivity and microbial degradation rates. Results of sensitivity studies show that including two ligands improves the model-data fit significantly compared to a model with only one ligand class or with constant ligand concentration. We use our model to investigate the role of organic ligands in controlling the fate of aerosol iron in a changing climate and future primary production; to isolate the effects we compare model runs with and without changes in ligand distribution.