Flows in variably saturated media that exhibit second-type heterogeneity, in which abrupt changes of medium parameters occur, are simulated by the Green element method (GEM). Such media are usually encountered where soil formations have arisen by different geological or geomorphological processes spread over different time scales. Two challenges are posed when simulating flows in multiply zoned unsaturated media: one is the highly nonlinear nature of the flow within each zone, and the other is dealing with sharp contrast in medium parameters at the interfaces of different zones. Both challenges are accommodated in this paper using a flux-based Green element formulation to simulate the flow and incorporating the Picard and Newton-Raphson (N-R) algorithms to simplify the nonlinear discrete equations. Calculations are carried out on three numerical examples of infiltration into unsaturated soils in two spatial dimensions. The convergence rate of the N-R algorithm is superior to the Picard algorithm only for the first example, while none of the algorithms has a clear advantage for the other two examples. The N-R algorithm suffers from repeated calculation of derivatives of the medium parameters with respect to the pressure head, thereby compromising the accuracy of the solution and increasing computational cost.