In general, the water content (θ) form of Richards’ equation is not used when modeling water flow through layered soil since θ is discontinuous across soil layers. Within the literature, there have been some examples of models developed for layered soils using the θ-form of Richards’ equation. However, these models usually rely on an approximation of the discontinuity at the soil layer interface. For the first time, we will develop an iterative scheme based on Newton’s method, to explicitly solve for θ at the interface between 2 soils within a numerical scheme. The numerical scheme used here is the Method of Lines (MoL); however, the principles of the iterative solution could be used in other numerical techniques. It will be shown that the iterative scheme is highly effective, converging within 1 to 2 iterations. To ensure the convergence behaviour holds, the numerical scheme will be tested on a fine-over-coarse and a coarse-over-fine soil with highly contrasting soil properties. For each case, the contrast between the soil types will be controlled artificially to extend and decrease the extent of the θ discontinuity. In addition, the numerical solution will be compared against a steady-state analytical solution and a numerical solution from the literature.