Abstract. A three-dimensional coupled physical–biogeochemical model was used to simulate and examine temporal and spatial variability of sea surface pCO₂ in the Gulf of Mexico (GoM). The model was driven by realistic atmospheric forcing, open boundary conditions from a data-assimilative global ocean circulation model, and observed freshwater and terrestrial nutrient and carbon input from major rivers. A 7-year model hindcast (2004–2010) was performed and validated against ship measurements. Model results revealed clear seasonality in surface pCO₂ and were used to estimate carbon budgets in the Gulf. Based on the average of model simulations, the GoM was a net CO₂ sink with a flux of 1.11 ± 0.84  ×  10¹² mol C yr⁻¹, which, together with the enormous fluvial inorganic carbon input, was comparable to the inorganic carbon export through the Loop Current. Two model sensitivity experiments were performed: one without biological sources and sinks and the other using river input from the 1904–1910 period as simulated by the Dynamic Land Ecosystem Model (DLEM). It was found that biological uptake was the primary driver making GoM an overall CO₂ sink and that the carbon flux in the northern GoM was very susceptible to changes in river forcing. Large uncertainties in model simulations warrant further process-based investigations.