A numerical model for the combined reaction, flow, mass and heat transfer close to gas bubbles in a continuous liquid phase is developed. The mathematical model consists of a set of highly nonlinear equations due to the strong coupling of the mass and heat transport and the reaction kinetics. Our simulation results predict that for a certain range of operating conditions, different steady states and mass transfer rates may be obtained, i.e., bubbles can be either ignited or extinguished. Furthermore, for a narrow range of conditions, three steady states are stable, i.e., "cold", "warm" and "hot" bubbles may co-exist. The first two steady states are due to the classical ignition/extinction behavior in analogy to cooled CSTRs, and the third state is caused by multiplicity of the enhancement factor. Bifurcation diagrams are constructed using an efficient pseudo-arc-length continuation method coupled with a powerful Newton-Krylov-subspace solver. Bifurcation maps are generated using only matrix-vector products, thus making the computation of these maps feasible (> 105 equations). Multiplicity of reacting rising bubbles has not been reported in the literature. Thus, the current study may explain several inconsistencies in experimental investigations, and it may serve as a basis for the rational design, optimization and control of bubble columns. (C) 2004 Elsevier Ltd. All rights reserved.