Massive blooms of Phaeocystis colonies usually occur in the Belgian coastal zone (BCZ) between spring and summer diatom blooms but their relative magnitude varies between years. In order to understand this interannual variability, we used the biogeochemical MIRO model to explore the link between diatom and Phaeocystis blooms and changing nutrient loads and meteorological conditions over the last decade. For this application, MIRO was implemented in a simplified 3-box representation of the domain between the Baie de Seine and the BCZ. MIRO was run over the 1989–2003 period using actual photosynthetic active radiation (PAR), seawater temperature and riverine nutrient loads as forcing. The water mass residence time was calculated for each box based on a monthly water budget estimated from 1993–2003 water flow simulations of the three-dimensional hydrodynamical model COHSNS-3D. Overall MIRO simulations compare fairly well with nutrient and phytoplankton data collected in the central BCZ but indicate the importance of the hydrodynamical resolution frame for correctly describing the extremely high nutrient concentrations and biomass observed in the BCZ. Analysis of model results suggests that while interannual variability in diatom biomass depends on both meteorological conditions (light and temperature) and nutrient loads, Phaeocystis blooms are mainly controlled by nutrients. Further sensitivity tests with varying N and P loads suggest that only N reduction will result in significantly decreased Phaeocystis blooms without negative affects on diatoms, while P reduction will negatively affect diatoms. Moreover, Atlantic nutrient loads play such a great role in BCZ enrichment that reduction of Scheldt nutrient loads only is not sufficient to significantly decrease phytoplankton blooms in the BCZ. It is concluded that future nutrient reduction policies aimed to decrease Phaeocystis blooms in the BCZ without impacting diatoms should target the decrease of N loads in both the Seine and the Scheldt rivers.