The Bay of Marseille (BoM), located in the northwestern Mediterranean Sea, is affected by various hydrodynamic processes (e.g., Rhône River intrusion and upwelling events) that result in a highly complex local carbonate system. In any complex environment, the use of models is advantageous since it allows us to identify the different environmental forcings, thereby facilitating a better understanding. By combining approaches from two biogeochemical ocean models and improving the formulation of total alkalinity, we develop a more realistic representation of the carbonate system variables at high temporal resolution, which enables us to study air–sea CO2 fluxes and seawater pCO2 variations more reliably. We apply this new formulation to two particular scenarios that are typical for the BoM: (i) summer upwelling and (ii) Rhône River intrusion events. In both scenarios, our model was able to correctly reproduce the observed patterns of pCO2 variability. Summer upwelling events are typically associated with a pCO2 decrease that mainly results from decreasing near-surface temperatures. Furthermore, Rhône River intrusion events are typically associated with a pCO2 decrease, although, in this case, the pCO2 decrease results from a decrease in salinity and an overall increase in total alkalinity. While we were able to correctly represent the daily range of air–sea CO2 fluxes, the present configuration of Eco3M_MIX-CarbOx does not allow us to correctly reproduce the annual cycle of air–sea CO2 fluxes observed in the area. This pattern directly impacts our estimates of the overall yearly air–sea CO2 flux as, even if the model clearly identifies the bay as a CO2 sink, its magnitude was underestimated, which may be an indication of the limitations inherent in dimensionless models for representing air–sea CO2 fluxes.