The hydrodynamics involved with the rise of air bubbles in shear‐thinning non‐Newtonian liquids in bubble columns was investigated using computational fluid dynamics (CFD). In bubble columns, the bubble size distribution (BSD) and gas holdup significantly affect the mass transfer rates and the reactor design. Consequently, the influence of superficial gas velocity, flow index, consistency index, sparger bubble size, and yield stress were analyzed on a local and global scale. The bubble sizes and the coalescence and breakage phenomena were incorporated using the homogeneous discrete method of the population balance model (PBM). The bubbles underwent coalescence during ascent and exhibited bimodal distribution. The radial homogeneity of the gas phase increased axially. Moreover, the zones of low liquid dynamic viscosity produced zones of a high gas holdup. Though the dual effect of viscosity on gas holdup was non‐existent, the gas holdup surged after 35.1 mPa s at the mid‐zone and remained constant thereafter. The noteworthy differences in simulation results put further emphasis upon the importance of sparger bubble size in CFD modelling. A decrease in the overall gas holdup and the axial air velocity with yield stress was observed. CFD simulation proved capable of providing results in reasonable agreement with experiments.