In this work we discuss an approach for uncertainty propagation through computationally expensive physics simulation codes. Our approach incorporates gradient information information to provide a higher quality surrogate with fewer simulation results compared with derivative-free approaches. We use this information in two ways: we fit a polynomial or Gaussian process model ("surrogate") of the system response. In a third approach we hybridize the techniques where a Gaussian process with polynomial mean is fit resulting in an improvement of both techniques. The surrogate coupled with input uncertainty information provides a complete uncertainty approach when the physics simulation code can be run at only a small number of times. We discuss various algorithmic choices such as polynomial basis and covariance kernel. We demonstrate our findings on synthetic functions as well as nuclear reactor models.