In this study, a mathematical model is developed for the cathode of PEM fuel cells, including multi-phase and multi-species transport and electrochemical reaction under the isothermal and steady-state conditions. The conservation equations for mass, momentum, species and charge are solved using the commercial software COMSOL Multiphysics. The catalyst layer is modeled as a finite domain and assumed to be composed of a uniform distribution of supported catalyst, liquid water, electrolyte and void space. The Stefan-Maxwell equation is used to model the multi-species diffusion in the gas diffusion and catalyst layers. Owing to the low relative species' velocity, Darcy's law is used to describe the transport of gas and liquid phases in the gas diffusion and catalyst layers. A serpentine flow field is considered to distribute the oxidant over the active cathode electrode surface, with pressure loss in the flow direction along the channel. The dependency of the capillary pressure on the saturation is modeled using the Leverette function and the Brooks and Corey relation. A parametric study is carried out to investigate the effects of pressure drop in the flow channel, permeability, inlet relative humidity and shoulder/ channel width ratio on the performance of the cell and the transport of liquid water. An inlet relative humidity of 90 and 80% leads to the highest performance in the cathode. Owing to liquid water evaporation, the relative humidity in the catalyst layer reaches 100% with an inlet relative humidity of 90 and 80%, resulting in a high electrolyte conductivity. The electrolyte conductivity plays a significant role in determining the overall performance up to a point. Further, the catalyst layer is found to be important in controlling the water concentration in the cell. The cross-flow phenomenon is shown to enhance the removal of liquid water from the cell. Moreover, a shoulder/channel width ratio of 1:2 is found to be an optimal ratio. A decrease in the shoulder/channel ratio results in an increase in performance and an increase in cross flow. Finally, the Leverette function leads to lower liquid water saturations in the backing and catalyst layers than the Brooks and Corey relation. The overall trend, however, is similar for both functions.