Surrogate modeling is computationally attractive for problems that require repetitive yet expensive simulations of PDEs, such as uncertainty propagation, deterministic design or inverse problems, where the main challenges are curse of dimensionality, data efficiency, uncertainty quantification and generalization, especially for problems with high dimensional input. To tackle those issues, we propose an image-to-image regression approach with deep dense convolutional encoder-decoder networks to learn accurate surrogate models for problems with up to 4225 input dimensions that is out of the reach of traditional surrogate models. A Stein variational inference method is scaled to modern Bayesian deep neural networks to further boost the regression performance and provide well-calibrated uncertainty estimate even with limited training data. We further explore how to incorporate the governing equations of the physical models into the loss/likelihood functions of the physics-constrained surrogates to completely avoid any simulation (or labeled data) while achieving similar accuracy with the date-driven surrogates. A conditional flow-based generative model with reverse KL-divergence loss without labeled data is trained to capture the predictive uncertainty. We also extend those methods to dynamic problems by treating time as extra input for dynamic multiphase flow. The efficacy of those methods are demonstrated in a variety of surrogate modeling and uncertainty quantification tasks in heterogeneous porous media flow.