We use subsurface characterization to describe porous media properties, such as permeability and porosity. However, the subsurface characterization is a demanding task. One of the challenges in the characterization is that we need to deal with a large dimension of the stochastic space. Typically, a dimensional reduction method, such as a Karhunen-Loeve (KL) expansion, is applied to the prior distribution in a Bayesian framework to make the characterization computationally tractable. Due to the large variability of the properties in the subsurface formations, it is worth localizing the sampling strategy. This will enable us to capture the local variability of rock properties more accurately. In this talk, we introduce the concept of multiscale sampling to localize the search in the stochastic space. In the Bayesian framework, we combine the new multiscale algorithm with a preconditioned Markov chain Monte Carlo (MCMC) algorithm. The new sampling algorithm decomposes the stochastic space in orthogonal complement subspaces, through a one-to-one mapping to a non-overlapping domain decomposition of the region of interest. The localized search is performed by a Gibbs sampler. In that the KL expansion is applied locally at the subdomain level. We demonstrate the effectiveness of the proposed sampling algorithm for the solution of an inverse elliptic problem. In a multi-MCMC study performed on Graphics Processing Units (GPUs), we show that the new algorithm clearly improves the convergence rate of the preconditioned MCMC method.