ABSTRACT The temporal and spectral evolution of gamma-ray burst (GRB) afterglows can be used to infer the density and density profile of the medium through which the shock is propagating. In long-duration (core-collapse) GRBs, the circumstellar medium (CSM) is expected to resemble a wind-blown bubble, with a termination shock, separating the stellar wind and the interstellar medium (ISM). A long standing problem is that flat density profiles, indicative of the ISM, are often found at lower radii than expected for a massive star progenitor. Furthermore, the presence of both wind-like environments at high radii and ISM-like environments at low radii remains a mystery. In this paper, we perform a ‘CSM population synthesis’ with long GRB progenitor stellar evolution models. Analytic results for the evolution of wind blown bubbles are adjusted through comparison with a grid of 2D hydrodynamical simulations. Predictions for the emission radii, ratio of ISM to wind-like environments, wind, and ISM densities are compared with the largest sample of afterglow derived parameters yet compiled, which we make available for the community. We find that high ISM densities of n ∼ 1000 cm−3 best reproduce observations. If long GRBs instead occur in typical ISM densities of n ∼ 1 cm−3, then the discrepancy between theory and observations is shown to persist at a population level. We discuss possible explanations for the origin of variety in long GRB afterglows, and for the overall trend of CSM modelling to over-predict the termination shock radius.