A mathematic model considering the interaction of bubbles was developed to simulate bubbles behavior in a continuous slab casting mold. The flow of molten steel is calculated with an Eulerian model coupled with a Lagrangian model for bubbles. Modified bubble coalescence and breakup models are established as user-defined functions of FLUENT software. Collisions between bubbles are described by means of hard-sphere model. The break-up mold is implemented with a break-up constrain based on the critical Weber number augmented a model for daughter bubble size distribution. The daughter bubble fraction is calculated with a tanh function, which has a U-shape distribution with lower bounds. Small bubbles were produced by the breakup of large bubbles. The predicted size and number of captured bubbles obtained with the mathematic model agree well with the plant measurements. The size and number distributions of small bubbles near the advancing solidification shell were obtained with simulation results. The results show that the bubble size and number at both the narrow and wide faces decrease with the increase of distance from the meniscus, and they get close to a constant in the deep region. Meanwhile, the initial bubble size has no obvious effect on the bubble distribution in the mold at the operation conditions in the present work. (C) 2019 Elsevier B.V. All rights reserved.