Accurately predicting polydisperse bubbly flow is a nontrivial task due to the complexity of the bubble number density function (NDF) and the strong dependence of the instantaneous bubble velocity on the bubble size and shape. To describe polydisperse bubbly flow, a joint mass-velocity NDF is adopted in this work. In the absence of mass transfer between phases and coalescence or breakage, the bubble mass is a conserved quantity from which the bubble size and shape can be found given the liquid pressure and surface tension. Quadrature-based moment methods (QBMM) are applied to solve numerically the kinetic equation of the joint NDF using the extended quadrature method of moments (EQMOM) coupled with an open-source incompressible Navier-Stokes solver for the liquid phase. Transport equations for the joint mass-velocity moments are derived from a kinetic equation for the joint NDF and closure is attained using a monokinetic NDF valid in the limit of small bubble Stokes number. The integer moments with respect to mass are used to reconstruct the continuous univariate NDF with EQMOM, while the joint mass-velocity moments are used to determine the bubble velocity as a continuous function of the bubble mass. The model is first applied to simulate a quasi-2-D bubble column with different aeration profiles and a narrow bubble size distribution in order to validate the approach with experimental data from the literature. Additional cases with a wide continuous bubble size distribution are used to show the ability of the modelling approach to describe polydisperse bubbly flows.