Reservoir history matching is a computationally expensive process, which requires multiple simulation runs. Therefore, there is a constant quest for more efficient sampling algorithms that can provide an ensemble of equally-good history matched models with a diverse range of predictions using fewer simulations. We introduce a novel stochastic Gaussian Process (GP) for assisted history matching where realizations are considered to be Gaussian random variables. The GP benefits from a small initial population and selects the next best possible samples by maximizing the expected improvement (EI). The maximization of EI function is computationally cheap and is performed by the Differential Evolution (DE) algorithm. The algorithm is successfully applied to a structurally complex faulted reservoir with 12 unknown parameters, 8 production and 4 injection wells. We show that the GP algorithm with EI maximization can significantly reduce the number of required simulations for history matching. The ensemble is then used to estimate the posterior distributions by performing the Markov chain Monte Carlo (McMC) using a cross-validated GP model. The hybrid workflow presents an efficient and computationally-cheap mechanism for history matching and uncertainty quantification of complex reservoir models.