Benthic oxygen fluxes consist of advective and diffusive terms. Both terms in the south-eastern North Sea exhibit a prominent annual cycle but with opposite variation patterns. To understand the driving mechanisms quantitatively, a novel 3-D benthic-pelagic coupled model resolving interactions among macrobenthos, bioturbation, oxygen consumption and carbon early diagenesis was applied to reconstruct the benthic states. Simulation results show a satisfactory agreement with field data and reveal that the benthic oxygen flux is determined by not only pelagic drivers but also by internal dynamics associated with the interaction between organic carbon and macrobenthos, and bedform morphodynamics. Variation of advective flux, characterized by summer-low and winter-high, is mainly driven by hydrodynamics and bedform morphodynamics, while variation of diffusive flux, featured by summer-high and winter-low, is a compound effect of pelagic and benthic drivers with a dominant control by macrobenthos through bioturbation. The role of bioturbation in benthic oxygen consumption is twofold: (i) on the one hand, it alters the particulate organic carbon (POC) distribution in surface sediments, thereby changing the availability of POC to oxygen consumption; (ii) on the other hand, it mixes oxygen down into sediments, thereby facilitating oxygen consumption. Our results indicate that the first role prevails in sandy seafloor characterized by energetic hydrodynamics, while the second role becomes increasingly important along with a weakening of bottom currents. We found that bioturbation-induced oxygen consumption contributes to more than 85% and 52% of the total benthic oxygen fluxes in muddy seabed and at a regional scale, respectively.