This thesis discusses and compares two approaches to solve parabolic partial differential equations with random input data. The stochastic problem is firstly transformed into a parametrized one by using finite dimensional noise assumption and the truncated Karhunen--Lo\`eve expansion. Both of these two approaches are applied onto this parametrized problem. The first approach considers the random coefficients to be their expected functions. The second approach, Monte Carlo discontinuous Galerkin (MCDG) method, randomly generates $M$ realizations of uncertain coefficients and obtains the numerical solution by averaging the $M$ solutions. These two approaches are compared in two numerical examples. The first example is a parabolic partial differential equation with random convection term in two dimensional case. The solution by the first approach moves in the direction determined by convection term, while the MCDG solution spreads towards the whole domain and diffuses faster than the solution by the first approach. The second example is a benchmark coupled flow and transport system in two cases. I first apply second-order kernel PCA to generate $M$ realizations of random permeability fields. They are used to obtain $M$ realizations of random velocity fields from the flow equation by DG. In the first approach, I take the convection term in the transport equation to be the average of the $M$ velocity realizations and solve the transport equation once by DG. The solution leaves the initial location completely in the direction of the average velocity field. In the second approach, I solve the transport equation for $M$ times corresponding to $M$ velocity realizations. The MCDG solution spreads toward the whole domain from the initial location and yields different speeds, spread and plume of the concentration of contaminant. Comparison shows that MCDG solution is more realistic, because it takes the uncertainty in velocity field into consideration. Besides, in order to correct overshoot and undershoot solutions caused by high level of oscillation in random convection realizations, I apply permeability realizations of lower resolution on finer meshes, and use slope limiter as well as lower and upper bound constraints. Additionally, I apply the $M$ realizations of permeability fields to a single-phase flow model on coarser meshes. The computing time is reduced significantly while accurate solutions are obtained. Finally, the future work is proposed.