Persistent organic pollutants (POPs) are of considerable concern due to their well-recognized toxicity and their potential to bioaccumulate and engage in long-range transport. This work describes the extension of the Community Multiscale Air Quality (CMAQ) modelling system to POPs with a focus on establishing an adaptable framework that accounts for gaseous chemistry, heterogeneous reactions, and gas-particle partitioning (GPP). The effect of GPP is assessed by implementing a set of independent parameterizations within the CMAQ aerosol module, including the Junge-Pankow (JP) adsorption model, the Harner-Bidleman (HB) organic matter (OM) absorption model, and the dual Dachs-Eisenreich (DE) black carbon (BC) adsorption and OM absorption model. Use of these descriptors in a modified version of CMAQ for benzo[a]pyrene (BaP) results in different fate and transport patterns as demonstrated by regional-scale simulations performed for a European domain during 2006. The dual DE model predicted 24.1% higher average domain concentrations compared to the HB model, which was in turn predicting 119.2% higher levels compared to the baseline JP model. It is found that the heterogeneous reaction of BaP with O-3 may decrease its atmospheric lifetime by 25.2% (domain and annual average) and near-ground concentrations by 18.8 %. Marginally better model performance was found for one of the six EMEP stations (Kosetice) when heterogeneous BaP reactivity was included. Further analysis shows that, for the rest of the EMEP locations, the model continues to underestimate BaP levels, an observation that can be attributed to low emission estimates for such remote areas. These findings suggest that, when modelling the fate and transport of organic pollutants on large spatio-temporal scales, the selection and parameterization of GPP can be as important as degradation (reactivity).