Numerical aspects of a pore scale model are investigated for the simulation of catalyst layers of polymer electrolyte membrane fuel cells. Coupled heat, mass and charged species transport together with reaction kinetics are taken into account using parallelized finite volume simulations for a range of nanostructured, computationally reconstructed catalyst layer samples. The effectiveness of implementing deflation as a second stage preconditioner generally improves convergence and results in better convergence behavior than more sophisticated first stage pre-conditioners. This behavior is attributed to the fact that the two stage preconditioner updates the preconditioning matrix at every GMRES restart, reducing the stalling effects that are commonly observed in restarted GMRES when a single stage preconditioner is used. In addition, the effectiveness of the deflation preconditioner is independent of the number of processors, whereas the localized block ILU preconditioner deteriorates in quality as the number of processors is increased. The total number of GMRES search directions required for convergence varies considerably depending on the preconditioner, but also depending on the catalyst layer microstructure, with low porosity microstructures requiring a smaller number of iterations. The improved model and numerical solution strategy should allow simulations for larger computational domains and improve the reliability of the predicted transport parameters. The preconditioning strategies presented in the paper are scalable and should prove effective for massively parallel simulations of other problems involving nonlinear equations.