We introduce a parallel rejection scheme to give a simple but reliable way to parallelize the Metropolis-Hastings algorithm. This method can be particularly useful when the target density is computationally expensive to evaluate and the acceptance rate of the Metropolis-Hastings is low. We apply the resulting method to quantify uncertainties of inverse problems, in which we aim to calibrate a challenging nonlinear geothermal reservoir model using real measurements from well tests. We demonstrate the parallelized method on various well-test scenarios. In some scenarios, the sample-based statistics obtained by our scheme shows clear advantages in providing robust model calibration and prediction compared with those obtained by nonlinear optimization methods.