Abstract. The isotopes of carbon (δ13C) and nitrogen (δ15N) are commonly used proxies for understanding the ocean. When used in tandem, they provide powerful insight into physical and biogeochemical processes. Here, we detail the implementation of δ13C and δ15N in the ocean component of an Earth system model. We evaluate our simulated δ13C and δ15N against contemporary measurements, place the model's performance alongside other isotope-enabled models and document the response of δ13C and δ15N to changes in ecosystem functioning. The model combines the Commonwealth Scientific and Industrial Research Organisation Mark 3L (CSIRO Mk3L) climate system model with the Carbon of the Ocean, Atmosphere and Land (COAL) biogeochemical model. The oceanic component of CSIRO Mk3L-COAL has a resolution of 1.6∘ latitude × 2.8∘ longitude and resolves multimillennial timescales, running at a rate of ∼400 years per day. We show that this coarse-resolution, computationally efficient model adequately reproduces water column and core-top δ13C and δ15N measurements, making it a useful tool for palaeoceanographic research. Changes to ecosystem function involve varying phytoplankton stoichiometry, varying CaCO3 production based on calcite saturation state and varying N2 fixation via iron limitation. We find that large changes in CaCO3 production have little effect on δ13C and δ15N, while changes in N2 fixation and phytoplankton stoichiometry have substantial and complex effects. Interpretations of palaeoceanographic records are therefore open to multiple lines of interpretation where multiple processes imprint on the isotopic signature, such as in the tropics, where denitrification, N2 fixation and nutrient utilisation influence δ15N. Hence, there is significant scope for isotope-enabled models to provide more robust interpretations of the proxy records.