Subglacial drainage models represent water flow at the ice--bed interface through coupled distributed and channelized systems to determine water pressure, discharge and drainage system geometry. While they are used to understand processes such as the relationship between surface melt and ice flow, the combination of the number of uncertain model parameters and their computational cost makes it difficult to adequately explore the high-dimensional parameter space and construct robust model predictions. Here, we develop Gaussian Process (GP) emulators that make fast predictions accompanied by uncertainty of subglacial drainage model outputs. Using a truncated principal component basis representation, we construct a GP emulator for daily representation of subglacial water pressure. We also explore emulation of scalar variables describing drainage efficiency and configuration. We train the emulators using ensembles of up to 512 simulations varying eight parameters of the Glacier Drainage System (GlaDS) model on a synthetic domain intended to represent an ice-sheet margin. The emulators make predictions ~1000 times faster than GlaDS simulations, with error <3% for the water pressure field and ~5--9% for drainage efficiency and configuration. We apply the emulators to explore the eight-dimensional input space by computing variance-based parameter-sensitivity indices, finding that three parameters (ice-flow coefficient, bed bump aspect ratio and the subglacial cavity system conductivity) explain 90% of the water pressure variance. The GP emulator approach described here is well-suited to integrate observational data with models to make calibrated, credible predictions of subglacial drainage.