New benchmarks are used to test two classes of discretization methods available in the literature to solve bivariate population balance equations (2-d PBEs), and the applicability of these mean-field equations to finite size systems. The new benchmarks, different from the extensions of their 1-d counterparts, relate to prediction of kinetics of mixing in particle phase under: (i) pure aggregation of particles, called aggregative mixing, and (ii) simultaneous breakup and coalescence of drops. The discretization methods for 2-d PBEs, derived from the widely used 1-d solution methods, are first classified into two classes. We choose one representative method from each class. The results show that the extensions based on minimum consistency of discretization perform quite well with respect to both the new and the old benchmarks, in comparison with the geometrical extensions of 1-d methods. We next revisit aggregative mixing using Monte-Carlo simulations. The simulations show that (i) the time variation of the extent of mixing in finite size systems has power law scaling with the system size, and (ii) the mean-field PBEs fail to capture the evolution of mixing for reduced population of particles at long times. The sum kernel limits the applicability of PBEs to substantially larger particle populations than that seen for the constant kernel. Interestingly, these populations are orders of magnitude larger than those at which the PBEs fail to capture the evolution of total particle population correctly. (C) 2010 Elsevier Ltd. All rights reserved.