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 CO$_2$ fluxes and seawater pCO$_2$ 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 pCO$_2$ 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 pCO$_2$ 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 CO$_2$ fluxes, the present configuration of Eco3M_MIX-CarbOx does not allow us to correctly reproduce the annual cycle of air–sea CO$_2$ fluxes observed in the area. This pattern directly impacts our estimates of the overall yearly air–sea CO$_2$ flux as, even if the model clearly identifies the bay as a CO$_2$ sink, its magnitude was underestimated, which may be an indication of the limitations inherent in dimensionless models for representing air–sea CO$_2$ fluxes.