Quantitative predictions of fluid flow and transport in porous media are often compromised by multi-scale heterogeneity and insufficient site characterization. These factors introduce uncertainty on the input and output of physical systems which are generally expressed as partial differential equations (PDEs). The characterization of this predictive uncertainty is typically done with forward propagation of input uncertainty as well as inverse modeling for the dynamic data integration. The main challenges of forward uncertainty propagation arise from the slow convergence of Monte Carlo Simulations (MCS) especially when the goal is to compute the probability distribution which is necessary for risk assessment and decision making under uncertainty. On the other hand, reliable inverse modeling is often hampered by the ill-posedness of the problem, thus the incorporation of geological constraints becomes increasingly important. In the thesis, four significant contributions are made to alleviate these outstanding issues on forward and inverse problems. First, the method of distributions for the steady-state flow problem is developed to yield a full probabilistic description of outputs via probability distribution function (PDF) or cumulative distribution (CDF). The derivation of deterministic equation for CDF relies on stochastic averaging techniques and self-consistent closure approximation which ensures the resulting CDF has the same mean and variance as those computed with moment equations or MCS. We conduct a series of numerical experiments dealing with steady-state two-dimensional flow driven by either a natural hydraulic head gradient or a pumping well. These experiments reveal that the proposed method remains accurate and robust for highly heterogeneous formations with the variance of log conductivity as large as five. For the same accuracy, it is also up to four orders of magnitude faster than MCS with a required degree of confidence. The second contribution of this work is the extension of the distribution-based method to account for uncertainty in the geologic makeup of a subsurface environment and non-stationary cases. Our CDF-RDD framework provides a probabilistic assessment of uncertainty in highly heterogeneous subsurface formations by combining the method of distributions and the random domain decomposition (RDD). Our numerical experiments reveal that the CDF-RDD remains accurate for two-dimensional flow in a porous material composed of two heterogeneous geo-facies, a setting in which the original distribution method fails. For the same accuracy, the CDF-RDD is an order of magnitude faster than MCS. Next, we develop a complete distribution-based method for the probabilistic forecast of two-phase flow in porous media. The CDF equation for travel time is derived within the efficient streamline-based framework to replace the MCS in the previous FROST method. For getting fast and stable results, we employ numerical techniques including pseudo-time integration, flux-limited scheme, and exponential grid spacing. Our CDF-FROST framework uses the results of the method of distributions for travel time as an input of FROST method. The proposed method provides a probability distribution of saturation without using any sampling-based methods. The numerical tests demonstrate that the CDF-FROST shows good accuracy in estimating the probability distributions of both saturation and travel time. For the same accuracy, it is about 5 and 10 times faster than the previous FROST method and naive MCS, respectively. Lastly, we propose a consensus equilibrium (CE) framework to reconstruct the realistic geological model by the inverse modeling of sparse dynamic data. The optimization-based inversion techniques are integrated with recent machine learning-based methods (e.g., variational auto-encoder and convolutional neural network) by the proposed CE algorithm to capture the complicated geological features. The numerical examples verify that the proposed method well preserves the geological realism, and it efficiently quantifies the uncertainty conditioned on dynamic information.