Fuel cells are electrochemical devices that directly convert the chemical energy of reaction of a fuel and an oxidant (usually hydrogen and oxygen) into electricity. Detailed measurements of properties such as the fluid flow distribution, species concentration, temperature distribution and local current densities are difficult to obtain during the operation of a fuel cell. Nonetheless, information like that is critical for improving the fuel cell performance, reducing cost and identifying possible failure mechanisms. Therefore, numerical models that are capable of producing this information would be very useful for design and optimization as well as improved fundamental understanding of transport processes in fuel cells. In this project a two-dimensional non-isothermal and non-isobaric computational model of a PEM fuel cell is formulated. The model is coupled with a computational fluid dynamics model for diffusive transport in the electrodes and convective transport in the gas flow channel, which is built principally upon the conservation laws for mass, momentum and energy. Then the following four phenomenological equations are also involved: the Stefan-Maxwell equation for the description of multi-species diffusion; the Bulter-Volmer equation for the description of first-order reaction kinetics; the Nernst-Planker equation for the description of proton flux through membrane; and the Schlögl equation for the description of water velocity in membrane. The focus of the project is placed on the study of partial hydration in the membrane. For this propose [sic], we develop our membrane model partly based on empirical relationships to account explicitly for water diffusion, pressure distribution and electro-osmotic drag; meanwhile, thermal effects are introduced implicity via the variation of membrane transport parameters as functions of temperature. Although the solution of the present model is beyond the scope of the project, we are trying to embed the general solution ideas in our modeling work by giving, for all computation domains, detailed boundary conditions and corresponding properties and parameters that are critical for the prediction accuracy of the model.