A comprehensive non-isothennal, 3D computational model for proton exchange membrane (PEM) fuel cells has been developed, and implemented into a computational fluid dynamic (CFD) code. The model allows parallel computing, thus making it practical to perform well-resolved simulations for large computational domains. The model accounts for convective and diffusive transport and allows prediction of the concentration of species. Distributed heat generation associated with the electrochemical reaction in the cathode and anode is included in the model. The model solves for the electric and ionic potentials in the electrodes and membrane, and the local activation overpotential distribution is resolved, rather than assumed uniform, making it possible to predict the local current density distribution more accurately. Maximum current densities are predicted under the land areas as a result of the dominant influence of ohmic losses over concentration losses on the activity at the catalyst layer. A parametric analysis shows that substantially different spatial distributions can be obtained by varying the asymmetry parameter with no noticeable change in the global current density and polarization curve. Changing the conductivity radically alters the current distribution by changing the relative influence of ohmic to activation overpotentials. (C) 2004 Elsevier B.V. All rights reserved.