Hydrogen embrittlement in martensitic press-hardened steel (PHS) presents a great challenge to the automotive industry that employs largely PHS in the critical parts of body-in-white. It is often a complex problem due to the coupling of martensitic transformation and hydrogen diffusion. Experimental studies often encounter difficulties in resolving the hydrogen diffusion in martensitic microstructure and its relation to hydrogen embrittlement in PHS. The present work develops a combined model coupling the phase field model for martensitic transformation together with a modified hydrogen diffusion model. We demonstrated, through both experiments and simulations, that hydrogen in PHS preferably resides in microstructure defects such as lath boundaries (LB) and prior-austenite grain boundaries (PAGB) where the hydrostatic stress and dislocation density remain high. Therefore, these boundaries are at risk from hydrogen embrittlement when the average hydrogen concentration of the material reaches a critical value. Furthermore, the ability of the model in resolving microscopic distribution allows us to predict fracture modes at any location by simply comparing the local concentration with the critical hydrogen concentration.