Subglacial topography is critically important for simulating ice-sheet evolution and projecting sea-level contributions. However, the subglacial topography of the Antarctic Ice Sheet is sparsely measured. Obtaining a gridded topography map used in ice-sheet simulations requires interpolating the measurements or inverting topography from observations of ice velocity and surface elevation. Traditional inverse methods based on the mass conservation law usually produce a single topography that is overly smooth and does not capture the non-uniqueness of the solutions to mass conservation. In this study, we develop a new method that combines geostatistical simulations with Markov chain Monte Carlo to stochastically generate different realizations of mass-conserving subglacial topography for regions with high ice velocity. This method uses a two-step approach that iteratively generates large-scale topography trends by random-field perturbations and simulates small-scale topography by geostatistical simulation. We test this method on Denman and Totten glaciers. The final topography ensemble shows significant elevation differences from BedMachine and presents large topographic uncertainty. This topography ensemble can be incorporated into ensemble ice-sheet modeling, allowing for the propagation of topographic uncertainty into the uncertainty in sea level projections.