Comment on hess-2021-15

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.

-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.
- -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 3 ... 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.
-  (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 2 .
-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?  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 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 timevariable, 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 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).