Review of Brandthorst et al. ‘Coupling saturated and unsaturated flow: comparing the iterative and the non-iterative approach‘

Major comments.

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.

Minor comments.

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 only)  in which you track what amounts of water go where during resizing.

l. 158, 221. What is point i? Do you mean point 1?

Figure 2. Perhaps use a different color for the horizontal line denoting the phreatic level. (Minor point, please don’t bother if this takes too much time.)

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’?

Figure 3. I was expecting a yes/no decision at the closure criterion, looping back to the start of the iteration in case of ‘no’.

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.

The manuscript presents two methodologies of coupling for dimensionnally heterogeneous modelling of subsurface flow in the unsaturated zone and in the saturated zone. The interest of such an approach is that important diminutions of computation times may be obtained compared to fully 3D modeling approaches. The main drawback is that accuracies of the simulations are dommaged, still in comparison with 3D modelling, and more or less along the considered cases and coupling methodologies. Given the very large computation times that may be encountered in fully mechanistic hydrological modeling at the wtareshed scale, this problem is of great interest for the community of hydrological modeling. The manuscript contains an important material in terms of numerical results and provides relevant hints to compare the two coupling stategies under concern. Nevertheless the presentation of the considered theories and numerical experiments lack of rigor, and the writting of the manuscript is not clear enough. In some places additionnal computations may even be needed. Thus I think it should be thoroughfully reworks prior to publication. I recommend to reject the paper in its present form, and to encourage the authors to resubmit after having completing and improving it.

General comments:

- The considered equations should be defined more rigourously and rewritten. For instance, the double time derivative term in equation (2) is a non-standard formulation of Richards equation (see for instance Gottardi and Venutelli, 1993). I guess that considerations related to the order on magnitudes of those two time derivatives may be used to justify the adopted formulation, but it should be explicited. Moreover, the use of the same notation S_y for the specific yield in equation (1), which has classically a clear and well identified physical meaning (drainage porosity), and for the iteratively computed, time variable fitting parameter used in the iterative method to handle recharge fluxes from the unsaturated zone is confusing. I think that the latter one should be expressed as the sum of the true specific yield S_y and a new additional term used for the purpose of the coupling between the saturated zone and the non-saturated zone. This would not imply new computation, but simply to rewrite some equations and rescale some results. I think that the added value of such more accurate notations in terms of clarity and of ease of physical interpretation would be important.

- Convergence studies for mesh refinement and time stepping strategy are not evocated as it should be the case in any study producing Computational Fluid Dynamics results. In some places it may imper the possibility to understand the comparative behaviours of the proposed test cases. For instance if we consider the comparison of accuracies of test 2 and test 3, in the present form of the manuscript it is impossible to say what comes from the differences of meshes and what comes from the different physics under concern (e.g.: homogeneous versus heterogeneous soil).

- To the knowledge of the reviewer, an important example of hydrological model that couples dimensionnaly heterogeneous descriptions of flow in the saturated zone and in the unsaturated zone is MIKE-SHE (e.g.: Graham and Butts, 2005), which is for instance included in recent international benchmarking efforts for physically based hydrological modeling (e.g.: Kollet et al., 2016). The fact that works related to MIKE-SHE do not appear in the references of the manuscript make me think that the bibliographical survey on which the presentation of the background of the study is done should be consolidated.

Specific comments:

- l 136-137: “a Neumann boundary representing net flux from precipitation and evapotranspiration” : with the source/sink term of the equation (2), it is possible to represent actual evapotranspiration distributed in time and space according to water avalaibility in the soil (see for instance Orgogozo et al., 2019) ; please discuss the limitation associated with an a priori estimation of the actual evapotranspiration directly embedded in the Neumann boundary.

- l 151 : “collect the computed recharge (i.e. flux leaving over the bottom boundary) and interpolate the 2D map of groundwater recharge.” : You mean collect all the computed regarges for all time steps of the 1D Richards model since the previous time step of coupling ? Should be clarified.

- l 157 : “Add (or subtract) a ratio r of this water to the recharge computed in the next time step.” you mean the next time step of coupling ? Should be explicited.

- l 159 – 171 : The proposed way of chosing the ratio r is difficult to accept. In case of water table elevation, the ratio r could be fitted to keep unchanged across the mesh resizing process the total amount of water contained in the part of the domain that stays unsaturated, while in case of water table lowering a ‘field capacity’ water saturation could be prescribed to the cells that experienced desaturation in order to compute a total water amount to be distributed in the new 1D mesh, with an associated proper r ratio? Here the formulae proposed for the computation of the ratio r seems somewhat arbitrary. For instance the point (1) “the groundwater table rise or fall is also effected by lateral flows” is already taken into account in the 2D groundwater model. Besides, “the unsaturated zone is really compacted by a rise of groundwater levels” does not sound physical at all.

- l 183 : “(ν) the iteration counter” : with which loop is related this iteration counter is unclear at this point (it could be for instance with the time stepping of the 1D Richards equation or with the coupling time steping)? Although it becomes clear afterward, it should be explicited here, at the first occurence of (ν).

- l 203 eq (10) (see the first general comment): According to the basic derivation of the diffusivity equation for unconfined aquifer, the specific yield is equal to the drainage porosity of the considered porous medium – although it seems that it might be different for more elaborated derivations, according to the literature cited by the authors. What is the physical interpretation of the variations the specific yield computed by eq (10)? Is there a theoritical reason why the iterating on the values of the specific yield field in the aquifer should lead to convergence? In case this is a purely empirical methodology, are there cases for which divergence may occure? Other questions : the value of the ‘physical’ S_y parameters that appears in the equation (1) is only the seed of the iterative process at the first time step of simulation, and do not appear directly anymore in the course of the simulation for the evaluation of S_y^v, right? That is what I understand from table 2 for instance. It should be clarified here.

- l 221 – 224 : “The source/sink terms q_lat,i have an effect on the recharge (‘R(Q lat = 0) ≠ R(Q lat ≠ 0)’), which due to the nonlinearity of Richards equation cannot be quantified.” However at step (4) (l 216 – 217), an updated R^ν is computed that takes into account the q_lat,i ^ν ? I don’t understand.

- l 229 – part 2.3 Activity score: Difficult to follow. Lack of explanations and of references. There is also a problem of structure: since ‘The parameters and the model output f are defined in Sections 3.4 and 4.4, respectively.’, this part 2.3 is not possible to understand by itself at this point of the reading. This part should be reworked so that the reader may understand why it is interesting to use the activity scores for the sensitivity study, and on what motivated the choices of the parameters x and the output f(x). By the way, in table 3 part 3.4, the variable K_GW and K_UZ seems not to be defined in the manuscript ? And why chosing S_y as a parameter of the sensitivity study while it is subjected to iterative evolution of its value along computation in the iterative method (l 218 – table 2)?

- l 267 : ‘3.1 Test case 1: 1D flow’ lack of a figure that presents the geometry, the boundary conditions and the meshes for each models.

- Figure 4 : Wrong title for y-axis (Precipitation-PET, not just PET)

- l 276 : Precise which 1D model (pure Richards I suppose ? )

- l 281 : “The groundwater domain is divided into a 2 × 2 grid. Each groundwater cell is assigned a 1D model.” Then the groundwater model is 2D with only 2 cells in each direction ? I don’t understand.

- l 285 : “Since there is no variability along the 8000 m side, flow is effectively 2D in this test case.” Then it is useless and misleading to present it as a 3D computation ; the figures and the discussions should be reshape for presenting directly the test case as a 2D one. The comparison of computation times is also questionnable : to deal with a 2D case in 3D increase tremendously (and artificially) the computation time with a fully mechanistic 3D model. Here some additionnal simulations (dealing with the 2D problem in 2D) are needed for making the comparison of computation times.

- l 290 : “[...] assigning a minimal initial pressure head of −1.25 m” ; you mean that -1.25 m is the pressure head at the top of the domain ? Please clarify.

- l 291 : “Monthly varying rainfall (Fig. 6) is used as Neumann boundary condition for the land surface”. More precision about these data would be useful – e.g.: are they synthetic ? Of which type of climate are they representative ?

- l 292 : Table 2 is not timely introduced ; since it contains information for the 3 test cases, it should be placed either in the beginning or at the end of the presentation of the considered test cases, but not at the middle.

- l 293 – 294 : “grid size ∆x = ∆y = 100 m and ∆z =0.1 m.”. It makes a form factor of 10³ ... Any convergence study done for the mesh refinement?

- l 296 : “With the flow problem being 2D this means that the entire domain is acutally covered by 1D models.” Nevertheless as far as I understood the proposed methodologies it would be exactly the same if the case was a 3D one? And I don’t understand to which extent a 1D approximation for a 2D problem would be essentially more “acute” than a 1D approximation for a 3D problem?

- l 306-307 : “three different soil units are distributed throughout the domain as depicted in Fig. 7.” More information is needed here. Is this distribution synthetic? How has it been acquired / built ? Of which type of soil (sand, loam, clay ...) each unit is representative ?

- l 308 : “averaged arithmetically” Any tests for the use of harmonic or geometric mean instead of arithmetic mean?

- l 309 : “ In the vertical direction a non-uniform grid is used with smaller grid sizes close to the surface and a total of 50 cells.” Please provide the minimum and maximum sizes.

- l 315 : “The 1D models are placed at the center of each zone.” How are laterally averaged the porous medium properties in each 1D models covering 10*8 cells laterally?

- l 327 : “The residual saturation S_r = θ_r /θ_s and the specific storage S_s are excluded from the analysis and set to 0.01 and 0.0015, respectively.” Why have they been excluded ? To be justified, or at least discussed.

- Table 3 : The parameters K_GW and K_UZ are appearing in the manuscript for the first time in this table. The notations used in table 3 and those used in the equations (especially (1) and (2)) should be the same, of at least explicitly related.

- l 329 -330 : “The horizontal spatial resolution is again ∆x = ∆y = 10 m, whereas the vertical resolution is ∆z = 0.1 m as in the 2D flow case.” Once again a convergence study must have been done to justify the use of this mesh with a form factor of 10².

- l 332 : “The time step sizes are the same as in the previously described test cases.” Any convergence study for justifying the use of the proposed time stepping policy?

- l 345 : “A visual comparison indicates that the coupling applied by Beegum et al. (2018) yields a comparable accuracy.” Why not plotting the results of Beegum in Figure 8?

- l 356 Table 4: This table contains information for all test cases and then it is not at the right place, being presented in a part specific to test case 1. Besides, since in test case 1 there is no lateral flux and thus no iteration in the iterative methods, I wonder why the iterative method has a wall time twice time more long than non-iterative method, while this later one include an addtionel step of remeshing? To be discussed.

- l 363-364 : “The results by Beegum et al. (2018) have a similar accuracy and shape as the results of the iterative coupling approach.” Why not plotting them in Figure 9?

- l 365 : “When considering the non-iterative model, it is notable that initial time steps are an issue [...]” Any convergence study on time step ? What happens if smaller time steps are used?

- l 366-368 : “Both of these issues may be related to the reference model essentially acting as a bucket without any plausible steady state solution (i.e. steady state for the groundwater model would have groundwater tables far above the top of the domain).” Then why not chosing lower values of precipitation , so that a steady state may be reached?

- l 374 : “All values [of S_y⁰] are smaller than the proposed value of 0.28, although the difference is less than 0.03.” How the proposed value of 0.28 has been choosen? Are there correlations between the S_y values and the state of the groundwaters (e.g.: water table altitude, lateral fluxes intensity)?

- l 379 : “Both coupling schemes show a good agreement with the fully integrated 3D model.” It is hard to understand why the matching between the fully 3D computation and the 2.5D ones is better here for this 3D heterogeneous test case than in the 2D homogeneous test case 2. I noted that in test case 3 a finer mesh is used than in test case 2. May be that convergence issues are at stake?

- l 389-391 : “Areas with larger differences appear at similar locations for both coupling schemes showing the largest deviations of up to ∆H GW = 1.5 m along the y = 800 m boundary”. Why are they such discrepancies, and why there? These points should be discussed here.

- l 394-396 : “Overall, the specific yield values are decreasing when the groundwater table is rising and increasing when the groundwater table is falling (roughly between t = 1100 d and t = 2200 d).” Once more, a careful discussion of the physical meaning of S_y and its variation is needed.

- l 414 : “Note that the specific yield in the iteratively coupled model is not the value used for the non-iteratively coupled model defined in Table 3 but the value calculated by the model during the simulation.” I don’t understand how it is possible to make a sensitivity analysis on a parameter that is not constant and specified prior to computation, but time-variable, calculated along computation?

- l 417 : “Acitivity”

- l 421-422 : “When looking at Eqs. 1 and 7, one sees that S_y can be eliminated which explains why there is no influence of S_y under these conditions.” You mean that dh/dt = 0 at extremas ? To be clarified.

- l 431-432 : “The average S_y value shows some smaller fluctuations, but overall it converges to a value around S_y = 0.17, which is a plausible value.”. This is a too short discussion of the value of this key parameter that controls the exchanges between the saturated zone and the unsaturated zone in the iterative method. It should be interpreted physically. It seems to potentially encompasses a non clearly identified list of physical phenomena.

- l 436-437 : “This means that the specific yield is mainly depending on the unsaturated zone parameters. This is reasonable as its intention is to represent the missing unsaturated zone in the groundwater model.” Somewhat strange. According to the basic derivation of the diffusivity equation for unconfined aquifers, the specific yield should be a property of the saturated zone (drainage porosity). So may be that if this parameters depends mainly on the properties of the unsaturated zone, it means that it is not, or not only, a specific yield (see the first general comment)?

- l 442 : “in the case of the iterative model even consistent.” I am not sure of what you want to say, please be more specific.

- l 448 : “On the contrary, using more models could help decreasing the discrepancies in the less accurate areas close to the no-flow boundary at y = 800 m which are most likely caused by the soil heterogeneities and the simplified recharge and specific yield pattern due to the zonation.” These discrepancies are important (~1,5m), and their causes must be carefully assessed. Additional numerical experiment with lower and stronger soil heterogeneities or various zonation startegies could help to ensure that the proposed diagnostic is correct. From my point of view stating that “As this is a general issue for these kind of models and does not relate to the presented coupling strategies themselves, we do not investigate it further.” is not sufficient, at least without any bibliographical references as it is at present.

- l 458-459 : “Therefore the results of the coupled model are on average more accurate even though this test case is more complex than the 2D flow case.” Meshes also are different, and without proper convergence studies the impact of this point may not be assessed. The convergence studies must be done, and used for consolidating the discussions.

- l 464 : “is constantly h_UZ = −1.25 m at ≥ 1.25 m above the groundwater table” This should be made clear sooner (see teh comment on l 290).

- l 467 : “Which parameter is dominating depends on the current flow conditions.” This should be discussed in more details.

- l 470 : “comformably” is not specific/quantitative enough.

- l 478-480 : “This is not the case in this model as we cannot calculate this effect properly and we therefore keep ∆H due to recharge fixed (see Eq. 7).” However in equations (7), (8) and (10), it is clear that there is an iterative procedure that involves ∆H_UZ^v and ∆H_GW^v that evolve at each iteration v? I don’t understand.

- l 481 : “In the end, the specific yield is not a physical quantity but a model parameter.”. This statement seems too general ; while it is clearly the case in the proposed modeling approach, it is not the case in all formulation of the diffusivity equation in unconfined aquifers. Overal all this paragraph should be rewritten to better discuss the meaning of the concepts that are specific to the proposed methodology with a wording that should not rise ambiguities between these concepts and previously existing concepts. For instance:

- l 485-486 : “The aim of the specific yield is to represent the missing unsaturated zone in the groundwater model” You are talking about what you called a specific yield in your model. I think that it should have another name that ‘specific yield’, this latter word designing a concept that do have physical meaning and that is related to the properties (drainage porosity) within the saturated zone in the basic form of the diffusivity equation for unconfined aquifers (see the first general comment).

This article proposes two approaches for coupling saturated and unsaturated zones in a hydrological model. An explicit representation of the 3D flows is time consuming and a 1D approach to represent the unsaturated zone with a 3D approach to represent the saturated zone is an alternative which however raise the question of the representation the interface between these two environments which evolves over time. This lead, among other feature to variation in the specific yield. Two innovative approaches, one iterative and the other non-iterative, are presented in this article and compared to a 3D reference model.

General comments :

Even if the article is based solely on synthetic data, the presentation of the approach is clear and the results are convincing both in terms of the quality of the results and the efficiency of the calculation times. This makes it an interesting article that deserves to be published. My main concern is the conclusions on the lack of sensitivity of the unsaturated zone. To my point of view the constant monthly inflow may lead to a steady recharge of the saturated zone. It might be possible that with pronounced short term drying-wetting cycles, the role of surface parameter might be stronger, which could be amplified by the coupling between surface flux and soil water content due to transpiration regulation. In Figure 14 we can see that soil parameter have their strongest effect when changing the inflow regime (except the model warming period). I think the conclusion could be tempered on the prominent sensitivity of KGW

specif comments :

The shape of the equation is unusual. Considering that θ=S(hp)*Φ, I have difficulty to understand how d(θ)/dt lead to left member of equation 2 (θ being the volumetric soil water content). Can the author give reference. The specific storage Ss is not clear (dS/dhp ?)

L154 : How the saturation determined into the new cells (is water mass in the unsaturated layer preserved?

L162 1 ratio and three terms. Not clear

In Table 3 : Are KGW and KUZ conductivity at saturation? Here there is a decoupling of K Values. What happen in area that might belonging to the two domains.

L393-396 Is the feature described here expected? May be can be addressed when discussing the yield values in part 4.5.

Figure 12: errors are located in particular areas. Is there some explanation (boundary conditions, but not everywhere, heterogeneity patterns)?

L397 not clear what the ration is

Section 4.4 Are the parameters leading to exfiltration (data removed ) cover particular domain, that might be meaningfull

Fig 14-16 : is t[a] correspond to year unit (you may consider t[y]