Comment on hess-2021-15

The contribution of the paper is purely computational. This can be seen from the test cases, all of which are highly artificial. I do not consider this a problem. The paper is focused, does not claim more than it delivers, and what it delivers is relevant for the HESS readership and substantial. The paper is well organized and generally clear. In the minor comments below I indicate where I was not sure if I understood everything.

The conclusion in the abstract that the more elaborate iterative coupling is better and more robust but slower than the non-iterative coupling is a bit underwhelming. Perhaps being more specific would make it more informative.
The overview in the Introduction of the various coupled unsaturated-saturated models, their coupling strategies and the role of the specific yield is insightful and thorough.
l. 124-131. If you assign the same 1D unsaturated zone columns to multiple grid cells of the underlying 2D groundwater model, I believe you are making the implicit assumption that the weather is uniform over the 2D extent of the groundwater level, or at least over the regions assigned to each of the 1D columns. If this is indeed the case, it may be good to mention this. l. 156-173. If I understand correctly, if the resizing of an unsaturated zone column shortens the column, you have to add the storage in the part of the column that is cut off to the amount of water in the saturated zone to conserve mass during the resizing operation. When you increase the size of an unsaturated column, some water is transferred from the saturated to the unsaturated zone. I believe this is done in steps 4 and 5. In that case, I believe the ratio r should normally be equal to 1 to ensure mass conservation. You give several reasons to deviate from this. Equation (6) gives the expression for r you used, capped at 1. If you consider the entire system, does capping r at a maximum of 1 whilst allowing values smaller than 1 not necessarily lead to mass losses when resizing all unsaturated columns?
Incidentally, for this reason I do not understand why Erdal et al. (2019) set r to 0 instead of 1, but I am not reviewing that paper here. But it seems to indicate that I missed something here, unless a shortening/lengthening of an unsaturated zone column led to a corresponding increase/decrease in its volumetric water content to keep the total amount of water present in the column unchanged.
Perhaps it is possible to have a 'before and after' figure, table, or water balance for the saturated and unsaturated domain (1 column  l. 203. Equation (10) forces the specific yield to conform to the simultaneous gain/loss of water and the resulting change in in the groundwater level. Essentially, the specific yield is no longer a model variable but the ratio of the calculated water gain per calculated groundwater level change. The water level change is calculated using the 1D Richards' equation, which does not have the specific yield as a parameter. You thereby to a certain extent impose the unsaturated zone model upon the groundwater model and adapt one groundwater model parameter (the specific yield) to match both models. Interesting approach.
l. 226. 'constant' or 'spatially uniform'? l. 271-272. You state that you have hydrostatic equilibrium as the initial condition (zero flow), but in the next sentence state that you state that in most of the profile, the matric potential is vertically uniform, which amounts to a profile with unit gradient flow.
l. 281-282: repeats l. 276-277. l. 325. Why did you assume uniform distributions for the van Genuchten parameters or their logarithms? Table 4. To increase the readability of the Table independently from the text I suggest to clarify that the calculation times are for a PC in the first two cases and a supercomputer in the third case.
Section 4.2 Could there be an effect of the way you calculate r (Eq. (6)) on the accuracy of the non-iterative 2D model? l. 394. The spread in values of the specific yield is well beyond the plausible range. But since the specific yield is not really a parameter in your approach I am not sure if that should be worrying or not.
l. 407-407. I recommend that you mention the fact that the model is not designed to handle overland flow when you describe the models in section 2. l. 461-465. You implemented unit gradient flow in the top of the profile and no flow (hydrostatic equilibrium) in the lower part of the profile, with an abrupt boundary between the two. In the first time step, this creates a rather hectic situation at the interface between these two regions, where Richards' equation needs to smoothen the transition and create some flow in the area with initially stagnant water. I can imagine this causes an initial jump in the importance of the van Genuchten parameters.
l. 481-487. In your approach, the specific yield is no longer even a model parameter. The discussion of its physically impossible values of the specific yield is OK, but the rest of the paragraph is a bit contrived.
l. 503-504. Unfortunately, the saturated hydraulic conductivity is also the parameter with the largest degree of spatial variation.