On the basis of Navier-Stockes equation and convection-diffusion equation, combined with surface tension and penetration models, the equations of moment and mass transfer between bubble and the ambient non-Newtonian liquid were established. The formation of a single bubble from a submersed nozzle of 1.0 mm diameter and the mass transfer from an artificially fixed bubble into the ambient liquid were simulated by the volume-of-fluid (VOF) method. Good agreement between simulation results and experimental data confirmed the validity of the numerical method. Furthermore, the concentration distribution around rising bubbles in shear thinning non-Newtonian fluid was simulated. When the process of a single ellipsoidal bubble with the bubble deformation rate below 2.0 rises, the concentration distribution is a single-tail in the bubble's wake, but it is fractal when the bubble deformation rate is greater than 2.0. For the overtaking of two in-line rising bubbles, the concentration distribution area between two bubbles broadens gradually and then coalescence occurs. The bifurcation of concentration distribution appears in the rear of the resultant bubble.