A population balance equation based dynamic gibbsite crystallization model, incorporating crystal growth, nucleation, and agglomeration, was solved using two different numerical techniques, namely, a fully implicit finite element method (FEM) and an explicit discretized population balance (DPB) numerical scheme. Unlike previous crystallization modeling approaches, the agglomeration model in the FEM framework was formulated based on the Safronov agglomeration equation and its partial differential equation (PDE) approximation [Bekker, A. V.; Livk, I. Agglomeration process modeling based on a PDE approximation of the Safronov agglomeration equation. Ind. Eng. Chem. Res. 2011, in press. The FEM numerical solution is implemented using a fully implicit Newton's method Galerkin finite element algorithm, which applies automatic Gear-type time-step and nonuniform adaptive mesh strategies to help optimal solution convergence. The dynamic FEM and DPB model predictions of isothermal gibbsite crystallization, using estimated kinetic parameters, are validated against experimental crystallization data that were generated under process conditions relevant to Bayer gibbsite crystallization. Additional modeling results obtained for a modified set of kinetic parameters, prolonged crystallization times, and more complex crystal size distributions show that the novel FEM-based crystallization modeling framework offers a more accurate and computationally efficient model solution than that based on the DPB approach.