The cooling of biscuits after baking can be considered as a random phenomenon. The Fourier field equation for transient heat transfer is used as the basis for the deterministic modelling of a single biscuit in a laboratory cooler. The Monte Carlo method is applied to the Fourier equation to predict the dispersion in biscuit surface temperature at the end of the cooling process. Theoretical solutions are also used to characterise the variability in biscuit temperature by considering biscuit cooling as a probabilistic thermal process. The mean and standard deviation in biscuit thermal properties, initial biscuit temperature and heat transfer coefficient are measured and are considered as the input random variables. The numerical output and theoretical solutions are compared with experimental results. Parameter studies are conducted with the probabilistic model to assess the influence of variability in thermal properties on final biscuit temperature variability and to suggest process improvement strategies. (C) 2002 Elsevier Science Ltd. All rights reserved.