Numerical simulation has been performed to investigate the bubbly flow in the liquid with vertical temperature gradient. The objective of this study is to establish an accurate numerical prediction program of gas-liquid two-phase flows in a vertical temperature gradient condition, whose mathematical model is formulated by the Eulerian-Lagrangian model. The present numerical results reveal the temperature mixing mechanism and the fluid dynamical characteristics induced by the bubbly flow in the liquid with stratified temperature. The effects of bubble radius, void fraction, and gas flow rate on bubble-driven convective flow are considered, too.