Free surfaces of large bubbles are involved in gas injection into liquid through a top-submerged lance. The Volume of Fluid (VOF) multiphase model is chosen to simulate these flows. However, the selection of an appropriate turbulence model remains a point of discussion. Therefore, three representative turbulence models, namely the Renormalization Group (RNG) k-epsilon model, the Reynolds Stress model (RSM) and the Large Eddy Simulation (LES) model, respectively coupled with VOF simulations, are studied. The results are compared with Particle Imaging Velocimetry (PIV) experiments. Afterwards, the selected turbulence model is used to study the two-phase flow characteristics. It has been found that the three turbulence models are able to capture the main averaged features of the flow such as the flow structures and average velocity field. However, the RNG and RSM models are not capable of estimating the fluctuating flow parameters such as the Reynolds stresses. In contrast, the LES model gives very good agreement with experimental data. The subsequent study on flow characteristics also showed that the 3D VOF-LES model can describe the bubble behavior and surface fluctuation. Resorting to Fast Fourier Transform (FFT) method, it has been proven that the bubble frequency is around 20 Hz and the turbulence in the bubbly region is between the single phase -5/3 and two-phase -10/3 power law in the inertial sub-range. Finally, the effect of the contact angle between the formed bubbles and outer lance wall on the flow has been estimated. The results showed that contact angle has limited influence on the flow. (C) 2018 Elsevier Ltd. All rights reserved.