Correctly describing oxygen reduction within the cathode catalyst layer (CL) in modelling proton exchange membrane fuel cell is an important issue remaining unresolved. In this paper we show how to derive an agglomerate model for calculating oxygen reactions by describing dissolved oxygen in the agglomerates using two independent random processes. The first one is the probability that an oxygen molecule, which dissolves in the ionomer film on the agglomerate surface, moves into and then remains in the agglomerates; the second one is the probability of the molecule being consumed in reactions. The first probability depends on CL structure and can be directly calculated; the second one is derived by assuming that the oxygen reduction is first-order kinetic. It is found that the distribution functions of the first process can be fitted to a generalised gamma distribution function, which enables us to derive an analytical agglomerate model. We also expend the model to include oxygen dissolution in the ionomer film, and apply it to simulate cathode electrodes. The results reveal that the resistance to oxygen diffusion in ionomer film and agglomerate in modern CL is minor, and that the main potential loss is due to oxygen dissolution in the ionomer film.