Because of the inherent multi-scale heterogeneity of porous media and the limitation of single-subject observed data, we propose to combine deterministic and stochastic techniques to simulate heterogeneity. We select a coastal plain sediment system as an example to demonstrate and verify this approach. Firstly, we apply transition probability matrix to determine and delineate the nonstationary unconformity, and combine hydro-stratigraphy analyses to establish the field/large-scale, deterministic stratigraphy model. Secondly, we apply fence diagrams and CPT data to infer the horizontal mean length of hydrofacies, and then build Markov chain models for each depositional system and simulate the local/intermediate-scale, stochastic hydrofacies model. Finally, we combine the stratigraphy and hydrofacies models to get a multi-scale heterogeneous model embedded with quantitative and qualitative observed data, with both deterministic and stochastic characteristics. In order to study the influence of uncertainty in model parameters on solute transport, we build multiple realizations of two types of heterogeneous model and use them to simulate groundwater flow and solute transport. The parameter sensitivity analysis shows the 1st and 2nd spatial moments of the contaminant plume increase with the lateral average length of hydrofacies.