It is proposed a new (hybrid) approach for the mathematical description of the black liquor burning process in an industrial recovery boiler furnace. The system is divided into four different regions and the concentration of each chemical substance in each region is calculated through direct minimization of the corresponding Gibbs free energy. The particulate formation is described separately through a neural network trained with industrial data. The resulting hybrid model reproduces satisfactorily black liquor burning data obtained from industrial sources, and also data reported in the open technical literature. (C) 2004 Elsevier Ltd. All rights reserved.