Bubble columns are complex multiphase reactors producing large contact areas for mass transfer and subsequent chemical reactions. For allowing a deeper insight into this complexity and providing strategies for process optimisation, CFD (computational fluid dynamics) has become a standard tool. However, accurate models are needed to reproduce correctly such dynamic systems. Numerical computations of reactive bubble columns are conducted based on a LES-Euler/Lagrange approach. For extending this approach beyond the point-mass assumption for the dispersed elements, a bubble oscillation model based on stochastic generation of eccentricity and motion angle is proposed. The effects of bubble dynamics are modelled in the interfacial forces accounting for instantaneous eccentricity. A dynamic Sherwood number is used to consider bubble dynamics in the mass transfer. The influence of chemical reactions in species consumption within the liquid phase was modelled through the enhancement factor. Large Eddy Simulation (LES) was used for calculating the fluid flow with momentum source terms for bubbles and modelling sub-grid-scale (SGS) turbulence in the continuous phase through a transport equation for the turbulent kinetic energy. The effect of sub-grid-scale (SGS) turbulence on bubble motion, and also SGS turbulence modification by bubbles (BIT) is considered, introducing an extra source term in a kSGS equation. This source term automatically accounts for turbulence dissipation as well as enhancement (i.e. BIT). Numerical simulations of the chemisorption process for CO2 bubbles rising in highly concentrated NaOH solutions were compared with data presented in the literature. A very good agreement was only found when applying the full bubble dynamics model.