A coupled biogeochemical-physical ocean model is used to study the seasonal and long–term variations of surface pCO2 in the North Atlantic Ocean. The model agrees well with recent underway pCO2 observations from the Surface Ocean CO2 Atlas (SOCAT) in various locations in the North Atlantic. Some of the distinct seasonal cycles observed in different parts of the North Atlantic are well reproduced by the model. In most regions except the subpolar domain, recent observed trends in pCO2 and air–sea carbon fluxes are also simulated by the model. Over the longer period between 1960–2008, the primary mode of surface pCO2 variability is dominated by the increasing trend associated with the invasion of anthropogenic CO2 into the ocean. We show that the spatial variability of this dominant increasing trend, to first order, can be explained by the surface ocean circulation and air–sea heat flux patterns. Regions with large surface mass transport and negative air–sea heat flux have the tendency to maintain lower surface pCO2. Regions of surface convergence and mean positive air–sea heat flux such as the subtropical gyre and the western subpolar gyre have a higher long–term surface pCO2 mean. The North Atlantic Oscillation (NAO) plays a major role in controlling the variability occurring at interannual to decadal time scales. The NAO predominantly influences surface pCO2 in the North Atlantic by changing the physical properties of the North Atlantic water masses, particularly by perturbing the temperature and dissolved inorganic carbon in the surface ocean. We show that present underway sea surface pCO2 observations are valuable for both calibrating the model, as well as for improving our understanding of the regionally heterogeneous variability of surface pCO2. In addition, they can be important for detecting any long term change in the regional carbon cycle due to ongoing climate change.