Rigorous prediction of species concentration distribution has been challenging in the simulation of bubble column reactors. In the present manuscript, a theoretical model of turbulent mass transfer is introduced to the numerical simulation of the process of chemisorption of CO2 into NaOH aqueous solution in a bubble column reactor. In this model the turbulent mass diffusivity was expressed as a function of the concentration variance (c(2))over-bar and its dissipation rate epsilon(c), which were modeled to close the turbulent mass transfer differential equation system. As a result, along with the fluid velocity distribution, the species concentration distribution in the bubble column can be rigorously predicted simultaneously. By using the present model, the experience relying assumption of a constant turbulent Schmidt, which has been usually adopted in a traditional species concentration simulation, can be avoided. An RNG k-epsilon method combined with a population balance model for predicting the bubble size distribution was used to predict the gas-liquid velocity in the turbulent flow. The simulated pH values and carbonate concentration were shown to be in satisfactory agreement with the experimental data in literature. The simulations showed that the turbulent mass diffusivity varied along axial and radial directions in the bubble column and the calculated turbulent Schmidt number was not a constant. The proposed approach can be especially useful when empirical turbulent mass diffusivity is not available for simulating a similar process in a bubble column reactor. (C) 2018 Elsevier Ltd. All rights reserved.