This paper studies numerically the convection of water vapour in snowpacks using an Eulerian-Eulerian two-phase approach. The convective water vapour transport in snow and its effects on snow density are often invoked to explain observed density profiles, e.g. of thin Arctic snow covers, but this process has never been numerically simulated and analysed in a systematic manner. Here, the impact of convection on the thermal and phase change regimes as a function of different snowpack depths, thermal boundary conditions and Rayleigh numbers is analysed. We find considerable impact of natural convection on the snow density distribution with a layer of significantly lower density at the bottom of the snowpack and a layer of higher density located higher in the snowpack or at the surface. We find that emergent heterogeneity in the snow porosity results in a feedback effect on the spatial organization of convection cells causing their horizontal displacement. Regions where the snowpack is most impacted by phase changes are found to be horizontally extended and vertically thin, `pancake'-like layers at the top and bottom of the snowpack. We further demonstrate that among the parameters important for natural convection, the snowpack depth has the strongest influence on the heat and mass transfer. Despite our simplifying assumptions, our study represents a significant improvement over the state of the art and a first step to accurately simulate snowpack dynamics in diverse regions of the cryosphere.