Coupling saturated and unsaturated flow: comparing the iterative and the non-iterative approach

Fully integrated three-dimensional (3D) physically based hydrologic models usually require many computational resources. For many applications, simplified models can be a cost-effective alternative. The 3D models of subsurface flow are often simplified by coupling a 2D groundwater model with multiple 1D models for the unsaturated zone. The crucial part of such models is the coupling between the two model compartments. In this work we compare two approaches for the coupling. One is iterative where the 1D unsaturated zone models go down to the impervious bottom of the aquifer, and the other one is non-iterative and uses a moving lower boundary for the unsaturated zone. In this context we also propose a new way of treating the specific yield, which plays a crucial role in linking the unsaturated and the groundwater model. Both models are applied to three test cases with increasing complexity and analyzed in terms of accuracy and speed compared to fully integrated model runs. The non-iterative approach is faster but does not yield a good accuracy for the model parameters in all applied test cases, whereas the iterative one gives good results in all cases. Which strategy is applied depends on the requirements: computational speed vs. model accuracy.


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 nonsaturated 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 : "(E) 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 (E).
-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 E is computed that takes into account the q lat,i E ? 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 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 0 ] 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 noniteratively 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 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).