We compare two approaches for coupling a 2D groundwater model with multiple 1D models for the unsaturated zone. One is non-iterative and very fast. The other one is iterative and involves a new way of treating the specific yield, which is crucial for obtaining a consistent solution in both model compartments. Tested on different scenarios, this new method turns out to be slower than the non-iterative approach but more accurate and still very efficient compared to fully integrated 3D model runs.