Abstract. Human exposure to mercury (Hg) is a cause of concern, due to the biomagnification of the neurotoxic species monomethylmercury (MMHg) in marine ecosystems. Previous research revealed that commercial fish species in the Mediterranean Sea ecosystems are particularly enriched in Hg, due to a combination of physical and ecological factors. Since the fate of Hg depends on the interactions among several biogeochemical and physical drivers, biogeochemical modeling is crucial to support the integration and interpretation of field data. Here, we develop and apply a coupled transport–biogeochemical–metal bioaccumulation numerical model (OGSTM–BFM–Hg) to simulate the biogeochemical cycling of the main Hg species (HgII, Hg0, MMHg, and DMHg) in seawater, organic detritus, and through the planktonic food web. The model is applied to a 3D domain of the Mediterranean Sea to investigate the spatial and temporal variability of methylmercury (MeHg) distribution and bioaccumulation and major uncertainties in Hg cycling. Model results reproduce the strong vertical and zonal gradients of MeHg concentrations related to primary production consistently with the observations and highlight the role of winter deep convection and summer water stratification in shaping the MeHg vertical distribution, including subsurface MeHg maximum. The modeled bioaccumulation dynamics in plankton food webs are characterized by a high spatial and temporal variability that is driven by plankton phenology and is consistent with available field data of Hg concentrations in plankton, as well as with other indicators, such as bioconcentration factors (BCFs) and trophic magnification factors (TMFs). Model results pointed out that the increment in water temperature linked to a decline of deep convection can cause an increase in water MeHg concentrations with cascading effects on plankton exposure and bioaccumulation.