Predicting mechanical stress contours on skin resulting from local tissue rearrangement surgeries is needed to design optimal treatment plans and avoid wound healing complications. Finite element (FE) simulations of skin tissues have been shown to be a reliable tool in preoperative planning, yet, a major obstacle in the creation of predictive software comes from the inherent uncertainty in material properties of biological materials, and the high computational cost of creating and calibrating virtual surgery models. In this study we build computationally inexpensive surrogates to easily predict stress profiles for arbitrary material parameters and a range of defect sizes in three reconstructive scenarios: advancement, transposition, and rotation flaps. The surrogates are built by first creating a training data set of FE simulations that cover the input space of experimentally-determined skin properties from the literature. A reduced order representation of the training data set is achieved via principal component analysis, and computationally efficient surrogates are then created through Gaussian Process (GP) regression. We show that the GP surrogates predict stress contours with relative errors that are on average 2% in the l2-norm compared to the high-fidelity FE models. We apply the GP surrogates to predict differences in the probability densities of stress contours between two different age groups undergoing the same procedure. By replacing nonlinear FE models with accurate yet inexpensive models that can be evaluated for any combination of human skin material parameters and a range of defect sizes, we aim to enable calibration and prediction of stress contours in individualized clinical cases in the near future.