We approach this problem by developing a physically realistic model of a hydrogen PEMFC, and link it to an efficient mathematical algorithm to perform system optimization. The model focuses on major transport processes in the gas channels and the membrane electrode assembly, and also takes into account various modes of water transport within the membrane, including diffusion, electro-osmotic drag, and hydraulic permeation. The resulting system of partial differential equations is fully discretized using a finite volume scheme leading to a large-scale non-linear system of equations which are solved simultaneously using an interior point optimization (IPOPT) solver. Parametric studies are performed to investigate effect of various design and operating conditions on the critical water management issues, and on overall system performance. In particular effects of inlet gas velocity, membrane thickness, inlet gas humidity, current density, flow arrangement and operating pressure are examined. Model results predict that thinner membranes result in more uniform water distribution and improved performance due to lower ionic resistance. Also, an internal water recirculation mechanism is observed within the system in the counter-flow arrangement, which gives better performance than the co-flow arrangement, especially for the cases of low pressure and low inlet gas humidity. Further, optimization of output power density with respect to the current density and the inlet gas velocity are performed. Current results indicate that our methodology of simultaneously solving the large-scale system of PEMFC model equations with IPOPT solver is a promising technique to examine complex cases for system optimization including multi-phase, non-isothermal transport processes and can be used for design under uncertainty for PEMFCs.