The density and flux of particles in the freeboard of a freely bubbling fluidized bed is a function of the size, aspect ratio and frequency of bubbles erupting at the surface of the dense bed. To predict the velocity and amount of material thrown into the freeboard, a previous model for the eruption of a single bubble (Yule and Glicksman, 1989, Proceedings Engineering Foundation International Fluidization Conference, Canada) was reformulated to deal with the eruption of a multiple bubbles undergoing coalescence at the bed surface. The solution includes the explicit mechanics of particle elements in the nose of the leading bubble as well as the particle elements in the wakes of bubbles trapped between coalescing bubbles. The motive force for the particle acceleration is the excess gas flow through the cavity at the bed surface. The gas flow, in turn is found from a numerical solution spanning the bed from the upper surface to the distributor. The solutions for a single erupting bubble, two vertically aligned bubbles, and three vertically aligned bubbles were obtained numerically. The predicted vertical mass flux of particles decreases exponentially with height. The smooth fall off of the predicted flux with height above the bed is due to the consideration of separate particle elements which have different initial conditions at different points around the nose of the bubble. The solution for the erupting single sphere agrees with the measurements of George and Grace (1978, A.I.Ch.E.Symp. Ser. 74). For a freely bubbling bed, the frequency of bubble eruptions and the probability of coalescence at the bed surface can be obtained from a statistical model of the bubbling process which has been used before to predict accurately the expansion of bubbling beds (Glicksman et al., 1991, Chem. Engng Sci. 46, 1561-1571). Combining the statistical model with the predictions for single eruptions results in a prediction of the flux and density distributions in the freeboard. Measured density distributions above a freely bubbling bed, containing horizontal tubes, are in general agreement with the model predictions. The distributions are found to be primarily controlled by the bubble size at the bed surface, the volume fraction of bubbles in the bed and the excess gas velocity U/U-mf. The predicted exponential coefficient for density and mass flux vs height is consistent with the correlation of experimental results given by Kunii and Levenspiel (1990, Powder Technol. 61, 193).