Abstract. We discuss several time integration techniques of system of ordinary differential equations characterized by strong nonlinear coupling of the unknown variables. This system is a result of spatial discretization (using such as finite element, finite difference, or finite volume element) of the Richards Equation which is a governing mathematical principle for modeling water infiltration through a subsurface. The nature of Richards Equation is further complicated by the fact that the rate of change of the quantity of interest represented by a time derivative is also nonlinear. We formulate a general framework of the numerical time integration as a discontinuous Galerkin method. The actual implementation of a particular scheme is realized by imposing certain finite element space in time variable to the variational equation and appropriate "variational crime" in the form of numerical quadrature for calculating the integration in the formulation. The resulting nonlinear algebraic equations are solved by employing some fixed point type iterations. We discuss two examples and compare their performance.