The thesis presents a new method for the solution of saturated-unsaturated groundwater flow problems in heterogeneous porous media. Concretely, highly nonlinear degenerate elliptic problems arising from a certain time discretization of the Richards equation are the basis of this work. The problems are considered as homogeneous in subdomains where a single soil prevails and, therefore, the parameter functions do not depend on space. These nonlinearities, however, may jump across the interfaces between the subdomains and, thus, account for the heterogeneous setting of different soil types in different subdomains. As a consequence, non-overlapping domain decomposition problems in which subproblems are coupled via nonlinear transmission conditions are obtained. In this work these problems are solved without any linearization. By Kirchhoff transformation the homogeneous subproblems are transformed into convex minimization problems. Here, additional constraints like Signorini-type boundary conditions, which occur on seapage faces around lakes, can be taken into account. Finite elements are chosen for the space discretization, and convex analysis is applied as the solution theory. Finally, monotone multigrid methods provide efficient solvers which are robust with respect to degenerating soil parameters. In order to deal with the coupling of the homogeneous subproblems, nonlinear Dirichlet-Neumann and Robin methods are used. Here, the thesis provides new convergence results for these iterations applied to nonlinear elliptic problems in 1D as well as well- posedness results, which generalize existing linear theory. On the other hand, detailed numerical experiments demonstrate that the methods can also be applied successfully to problems in 2D. Finally, based on the artificial viscosity method, an upwind discretization with finite elements is developed in order to account for gravity. Hence, stability of the numerical solutions is obtained. In a closing numerical example the Richards equation is solved in 2D with four different soils and coupled to a surface water reservoir. The result demonstrates the applicability of the developed solution technique to a heterogeneous problem with realistic hydrological data.