In this work we use computational fluid dynamics (CFD) to simulate a reactive gas-liquid bubbly system in a rectangular bubble column, operating at low superficial velocities (i.e. homogeneous regime). The gas bubbles, injected in the column through a sparger, contain one of the reactants, namely CO2, that via mass transfer moves to the continuous liquid phase, where it reacts with NaOH. A key role is played by the bubble size distribution (BSD) and the specific surface area that define the overall mass transfer rate in the CFD model. In order to correctly predict the BSD and the polydispersity of the bubbly system the population balance equation is solved by the quadrature method of moments (QMOM), within the OpenFOAM (v. 2.2.x) two-fluid solver compressibleTwoPhaseEulerFoam. To reduce the computational time and increase stability, a second-order operator-splitting technique for the solution of the chemically reactive species is also implemented, allowing to solve the different processes involved with their own time-scale. To our knowledge this is the first time that QMOM is employed for the simulation of a real reactive bubbly system and predictions are validated against experiments