Fuel cells are electrochemical devices that rely on the transport of reactants (oxygen and hydrogen) and products (water and heat). These transport processes are coupled with electrochemistry and further complicated by phase change, porous media (gas diffusion electrodes) and a complex geometry. This thesis presents a three dimensional, non-isothermal computational model of a proton exchange membrane fuel cell (PEMFC). The model was developed to improve fundamental understanding of transport phenomena in PEMFCs and to investigate the impact of various operation parameters on performance. The model, which was implemented into a Computational Fluid Dynamics code, accounts for all major transport phenomena, including: water and proton transport through the membrane; electrochemical reaction; transport of electrons; transport and phase change of water in the gas diffusion electrodes; temperature variation; diffusion of multi-component gas mixtures in the electrodes; pressure gradients; multi-component convective heat and mass transport in the gas flow channels. Simulations employing the single-phase version of the model are performed for a straight channel section of a complete cell including the anode and cathode flow channels. Base case simulations are presented and analyzed with a focus on the physical insight, and fundamental understanding afforded by the availability of detailed distributions of reactant concentrations, current densities, temperature and water fluxes. The results are consistent with available experimental observations and show that significant temperature gradients exist within the cell, with temperature differences of several degrees Kelvin within the membrane-electrode-assembly. The three-dimensional nature of the transport processes is particularly pronounced under the collector plates land area, and has a major impact on the current distribution and predicted limiting current density. A parametric study with the single-phase computational model is also presented to investigate the effect of various operating, geometric and material parameters, including temperature, pressure, stoichiometric flow ratio, porosity and thickness of the gas diffusion layers, and the ratio between the channel with and the land area. The two-phase version of the computational model is used for a domain including a cooling channel adjacent to the cell. Simulations are performed over a range of current densities. The analysis reveals a complex interplay between several competing phase change mechanisms in the gas diffusion electrodes. Results show that the liquid water saturation is below 0.1 inside both anode and cathode gas diffusion layers. For the anode side, saturation increases with increasing current density, whereas at the cathode side saturation reaches a maximum at an intermediate current density (≈ 1.1Amp/cm2) and decreases thereafter. The simulation show that a variety of flow regimes for liquid water and vapour are present at different locations in the cell, and these depend further on current density. The PEMFC model presented in this thesis has a number of novel features that enhance the physical realism of the simulations and provide insight, particularly in heat and water management. The model should serve as a good foundation for future development of a computationally based design and optimization method.