A frequent problem during numerical computations consists in the uncertainty of certain model parameters due to measuring errors or their high variability. In the last years, one could observe an increasing interest in the quantification of these uncertainties and their effects to the solution of numerical simulations; a powerful tool which has been proven to be an efficient approach in this context is the so-called polynomial chaos method which is based on a spectral decomposition of the covariance function of the uncertain parameters and a representation of the solution in a polynomial basis. The aim of this thesis is the application of this method to the Richards equation modeling groundwater flow in saturated and unsaturated porous media. The main difficulty consists in the saturation and the hydraulic conductivity appearing in the time derivative and in the spatial derivatives, since both depend nonlinearly on the pressure. Considering uncertain parameters like random initial and boundary conditions and, in particular, a random permeability leads to a stochastic variational inequality of second kind with obstacle conditions and a nonlinear convex functional as superposition operator. Considering variational inequalities in the context of uncertain parameters and the polynomial chaos method is new, and we start by deriving a weak formulation of the problem and approximating the parameters by a Karhunen-Loève expansion. The existence of a unique solution u in a tensor space can be proven for the time-discrete problem by reformulation as a convex minimization problem. We proceed by discretizing with finite elements and polynomial ansatz functions and by approximating the convex functional with Gaussian quadrature. The convergence of the solution of the discretized problem to the solution u is proved in a special case for a stochastic obstacle problem. Moreover, we perform numerical experiments to determine the discretization error. In the second part of this thesis, we develop an efficient numerical method to solve the discretized minimization problems. It is based on a global converging Block Gauß Seidel method and exploits a transformation which decouples the stochastic coefficients and connects the stochastic Galerkin with the stochastic collocation approach. This also allows us to establish a multigrid solver to accelerate the convergence. We conclude this thesis by demonstrating the power of our approach on a realistic example with lognormal permeability and exponential covariance.