In this thesis, two numerical modeling methods are used to investigate the thermal conductivity of the polymer electrolyte membrane (PEM) fuel cell gas diffusion layer (GDL). First, an analytical model is used to study the through-plane thermal conductivity from representative physical GDL models informed by microscale computed tomography imaging of four commercially available GDL materials. The effect of the heterogeneity of the through-plane porosity of the GDL and polytetrafluoroethylene (PTFE) treatment is studied and it is noted that the high porosity surface transition regions have a dominating effect over the addition of PTFE in impacting the overall thermal conductivity. Next, the lattice Boltzmann method (LBM) is employed to study both the in-plane and through-plane thermal conductivity of stochastic numerically generated GDL modeling domains. The effect of GDL compression, binder content, PTFE treatment, addition of a microporous layer (MPL), heterogeneous porosity distributions, and water saturation on the thermal conductivity are investigated.