Unlike the traditional two-stage methods, a conditional and inverse-conditional simulation approach may directly generate independent, identically distributed (i.i.d) realizations to honor both static data and state data in one step. The Markov chain Monte Carlo (McMC) method was proved a powerful tool to perform such type of stochastic simulation. One of the main advantages of the McMC over the traditional sensitivity-based optimization methods to inverse problems is its power, flexibility and well-posedness in incorporating observation data from different sources. In this work, an improved version of the McMC method is presented to perform the stochastic simulation of reservoirs and aquifers in the framework of multi-Gaussian geostatistics. First, a blocking scheme is proposed to overcome the limitations of the classic single-component Metropolis-Hastings-type McMC. One of the main characteristics of the blocking McMC (BMcMC) scheme is that, depending on the inconsistence between the prior model and the reality, it can preserve the prior spatial structure and statistics as users specified. At the same time, it improves the mixing of the Markov chain and hence enhances the computational efficiency of the McMC. Furthermore, the exploration ability and the mixing speed of McMC are efficiently improved by coupling the multiscale proposals, i.e., the coupled multiscale McMC method. In order to make the BMcMC method capable of dealing with the high-dimensional cases, a multi-scale scheme is introduced to accelerate the computation of the likelihood which greatly improves the computational efficiency of the McMC due to the fact that most of the computational efforts are spent on the forward simulations. To this end, a flexible-grid full-tensor finite-difference simulator, which is widely compatible with the outputs from various upscaling subroutines, is developed to solve the flow equations and a constant-displacement random-walk particle-tracking method, which enhances the computational efficiency at different scales, is employed to solve the transport problems. Second, the usefulness and efficiency of the proposed method are validated by a synthetic example. The uncertainty reduction due to conditioning on various types of data from different sources is assessed with the aid of the synthetic example. One of the novel achievements in this work is that the physical models are constrained to the temporal moments of BTCs that are more easily accessible than the concentration data which are only sparsely distributed in space. The worth on uncertainty reduction is evaluated by comparing to other data sources. Third, by comparing the BMcMC to the ensemble Kalman filtering (EnKF), the importance of honoring the prior information for inverse stochastic modeling is ascertained in two synthetic examples. Numerical simulations show that, even though the EnKF method may efficiently provide a better reproduction of observed dynamic data than the BMcMC method, the preservation of spatial statistics and model structure makes the BMcMC simulations competitive for some cases in predicting accurately and reliably the future performance of reservoirs particularly at new well locations. This is because the spatial structure and statistics of models may be one of the most important error sources to the prediction of the future performance of reservoirs and aquifers, it should be consistent with the given information just as conditioning to linear data and inverse-conditioning to nonlinear data. In other words, the realizations generated should preserve the given spatial structure and statistics during the procedure of conditioning and inverse-conditioning.