Numerical reservoir simulations are the most robust tool available for planning the sustainable utilization of geothermal resources. To increase the accuracy of simulation results, uncertainties in the input parameters and the essential driving factors of the power potential in fields under development need to be quantified via stochastic approaches and specified by sensitivity analyses, respectively. These quantifications and specifications are made here by applying a response surface methodology framework to a numerical dual-porosity model to generate probability distributions of the power potential and production decline rate and by analyzing pairwise parameter relationships. The Bedugul geothermal field (BGF) in Bali, Indonesia, with a complex structure composed of a vapor-dominated zone formed above a deep liquid reservoir, is chosen as a model case. The parameters selected as essential are permeability, porosity, water saturation, feed zone, volume factor, and reinjection because they control the power output and available period of geothermal resources. Using a multivariate regression equation, the uncertainties of six essential parameters are assessed using 54 proxy models built using the Box–Behnken design and the probability and cumulative distributions of Monte Carlo simulations of the models. The BGF is probabilistically revealed to contain a power potential of 92 MWe for 30 years of utilization, and the most influential parameters are the permeability for the reservoir temperature at the natural state and the feed zone depth for the power potential and decline rate.