Review of: Controls on managed aquifer recharge through a heterogeneous vadose zone: hydrologic modeling at a site characterized with surface geophysics, by Perzan et al

Summary and Recommendation

The paper describes a methodology for inferring coarse vs. fine grained fractions of the subsurface based on towed subsurface resistivity mapping (tTEM) calibrated to cone penetration test (CPT) data, in an almond orchard in the Central Valley California. Using various spatial and multi-variable statistics, this mapping is transformed to subsurface coarse-grained-fraction 3D realizations. Variably saturated flow model that simulates managed aquifer recharge by flooding part of the orchard in the spring is used to test hydrological (infiltration and recharge) and agricultural (root-zone saturation periods) farmers' interest. Results showed that coarse-grained structures accommodate the rapid recharge. Fine-grained sediments in the root zone are probable to cause long –term saturation affecting yield. Fine grained-blocks in the deeper unsaturated zone may retain significant volumes of the MAR operation water for years after the flooding. Infiltration, recharge and root-zone saturation were found sensitive only to the fine-grained fraction's hydraulic properties, therefore the authors conclude that better characterization of the fine-grain end member is needed to reduce uncertainty in similar Agricultural-MAR operations.

The innovative technical procedure including non-invasive and invasive geophysical methods, geostatistical, Monte-Carlo and multivariate statistics, and 3D variably saturated flow simulations, is of new and advanced nature and very well described, therefore a good fit for publication in HESS. I have some arguments on the overall understanding of recharge under thick unsaturated zone and the practical conclusion for Ag. MAR operation which I would like the authors to discuss and perhaps rethink. Therefore, I recommend moderate revisions (betweem minor and major).

Major Comments

1)    As the authors describe nicely recharge is correlated good with surface input (precipitation, irrigation, flood-MAR) due to the pressure wave that propagates fast in the unsaturated zone and push deep-unsaturated older water to the water-table (sometimes called by old groundwater hydrologists as "the train-car model"). Today we can deal with the unsaturated-zone's retention and avoid the simplistic recharge coefficients used in many groundwater models. Nevertheless, in the life-time of an almond-grove (20-50 years) the result of retention of ~ 1/3 of the MAR water in the deep unsaturated zone after 2 years, is no news and hardly important and shouldn’t be highlighted as a main result in the abstract. Higher storage in the deep unsaturated-zone will increase recharge in 1 of the following years of high input either due to a rainy year or MAR operation. The authors should discus and maybe reconsider the important implication of their analysis in the scope of Ag–MAR in orchards.  

More-ever, recent 3D simulations in heterogeneous variably saturated medium concerning transient flow from a drywell showed that fine-grained layers at the bottom of a dry well contribute to faster downward flow in the unsaturated zone and faster recharge - see Russo et al., 2022, WRR, https://doi.org/10.1029/2021WR031881. This phenomena is due to the turn over in unsaturated hydraulic conductivity during drying, where at increasing negative pressure head, fine-grain medium becomes more permeable than course-grain.

2)    Although used often, the term "deep vadose-zone" is awkward, as vadose comes from Latin meaning shallow. I suggest to use deep unsaturated zone for the domain between the bottom of root zone and the water table (especially if it is of tens of meters thick).

Specific comments

1)    L71 see also Rudnik et al., 2022, WRR for use of stochastic approaches in MAR

2)    L76-77 as discussed in major comment 1, in transient heterogeneous unsaturated flow fine-grained layers can increase flow in drying periods. Perhaps change "restrict flow" to impact flow.

3)    L113-126,Why not define the recharge of a 40 m deep aquifer as the downward flux at 39 m depth and stay with fluxes:1) it is straightforward and simpler; 2) the saturated zone part of a model may include sources and sinks (pumping wells) or transient head boundary conditions which have impact on the lateral flow not related to the MAR operation. Discuss.

4)    Figure 3, caption last row, change "hydraulic conductivity" to saturated hydraulic conductivity  

5)    L 209-210, perhaps better: Algebraically the resistivity of a tTEM cell is described by the harmonic mean of the fine and coarse grain resistivity's as:

6)    L249, May be better hydraulic functions than "water retention curve" (to include the unsaturated hydraulic conductivity function as well as the retention curve).

7)    L 360—370, only complete saturation? Or defined from some threshold of high saturation (e.g. 95% saturation)?

8)    L414 – 415 "or from the particular …simulation."   not clear, explain or discard if not important for the sensitivity analysis description.

9)    L 438 "0.15 +- 0.29" 0.29 standard deviation? define explicitly

10) L 461 discard "( 7 acres)"

11) Figure 6 caption. Is it a single flooding of 0.8 m, or another scheme? Should be said in caption.

12) Figure 7: 1) What drives the Flux recharge after 6 years when the water table is back to its pre-MAR level, perhaps not a consequence of MAR (relates to specific comment # 3); 2) right hand vertical axis title - typo recharge.

13) L494-495, this is a trivial result no need for numbers and statistics.

14) L 540 "especially within the vadose zone", where else than the unsaturated zone?

15) L616-621 – conclusion 2 – As said in major comments, the significance of this result in Ag.-MAR in almond groves is minor. 63% in 2 years for free or cheap water is no good? And the rest 37% are not lost forever they will recharge in the next rainy/MAR year (unless the pre-MAR unsaturated zone was in really low water contents)

16) L 629-634 conclusion 5 – Hard to belief that 1 flooding of 80 cm did not cause more saturation in rootzone than 16 inundations of 5 cm with 1 week between inundations. Check! A weekly 5 cm irrigation should not cause saturation of the entire root-zone unless very poor percolation in the soil.   

17) L633 saturated hydraulic conductivity rather than "hydraulic conductivity"     

Summary

Perzan et al. manuscript deals with modeling managed aquifer recharge (MAR) under an almond orchard using geophysical data and Monte-Carlo simulations to deal with a relatively deep heterogeneous vadose zone. The manuscript is well written, the figures are excellent, and the scope is well suited for HESS. The workflow of parametrizing the model using geophysical data is thoroughly done, and I like the authors’ approach to quantifying MAR recharge quantity and dynamics (comparing simulations with and without MAR). Unfortunately, the authors did not provide any validation of the simulation results by measurement of groundwater level and water content.

To give a wide perspective of their work, I think the authors tried to present this case study as representing a typical flood-MAR operation, that can include both agricultural MAR (Ag-MAR) and dedicated MAR facilities. This is not the case, in my opinion, and I suggest that the authors change the focus of the paper (including the title) to “Ag-MAR” or if they prefer “on-farm flooding”. This is not just semantics: dedicated MAR facilities usually have a limited spreading area that overlies high conductivity subsurface layers, while in Ag-MAR the tradeoff is an “unlimited” spreading area that overlies less conductive sub-surface. The submitted manuscript exactly deals with the second case.

I have some concerns with the modeling of the infiltration process, calculation of the root zone residence time (discretization and root zone sizes are similar), and some suggestions for improving the manuscript (see details below). Overall, I recommend minor revisions.

General comments

I strongly suggest revising the Methods section. It is too long - about half of the manuscript. Although very informative, some parts of it should be moved to the supplementary material, mainly parts that were published previously (e.g., Fig. 4 and related text which is a similar approach to Goebel and Knight, 2021).

Please elaborate more on the impact of model discretization on the infiltration process. The vertical discretization of your model is very rough (1 m) for numerical infiltration problems where usually discretization is on the order of 0.05-0.3 m (e.g., Botros et al. 2012). Moreover, this may have a profound impact on your root zone residence time analysis – because your vertical discretization (1 m or more) is the same size as the root zone (1 m). Hence, the root zone in your model is probably only one cell that has no spatial dynamics – at each time step it always has only one value of water content. Please explain if this can impact your analysis (e.g., overestimating root zone residence time).

Specific comments

L76-77: This is the motivation of the study, but it is not convincing, as other studies explore “The influence of meter-scale heterogeneity on recharge..” – I suggest to mention some of these studies and their findings.

L95: The authors chose the term flood-MAR to describe their framework, but it may be too general for their specific agricultural flood-MAR operation. Their choice of 2^nd metric of root zone residence time demonstrates this point. I suggest using the term “Ag-MAR” (e.g., Ganot & Dahlke, 2021; Levintal et al. 2022) or “on-farm flooding” (Bachand et al. 2014) instead of flood-MAR for most places in the text.

L111-112: In SAT-MAR anoxic conditions are one of the ways to reduce nitrate loads by denitrification.

L135-137: It is unclear why the authors chose for this work an agricultural field with a SAGBI rating of “poor” (maybe because of the available geophysical data?). This is an agricultural field that is not suitable for Ag-MAR according to the approach of O’Geen et al. (2015) which is much less exhaustive than the approach presented in the current paper. Please explain.

Table 1: Check the residual saturation (# 11 and 12). Values are too high (above common porosity values) – I believe it should be relative residual saturation.

L319-321: Not sure how old is the almond orchard at your site, but during very long spin-up times (>15-25 years) other crops grown at the site with different Kc may have a different impact on water budget and subsurface storage. In other words, it seems like very long spin-up times (such as 131 years) have a little practical benefit.

L335-342: The falling head boundary is elegant, but probably far from real-life water application of a large plot with an area of 800 m x 400 m. This is mainly true for an impractical initial condition with a ponding head of 0.8 m (!). From my experience ponding during on-farm flooding rarely exceeds 10 cm. Unrealistic high ponding conditions may also overestimate infiltration rates. An alternative upper boundary will be to keep a constant head of 5 to 10 cm followed by a falling head (e.g., for a total application of 20 cm, with a 5 cm constant head, let 15 cm infiltrate under a constant head, and then start the falling head loop for the remaining 5 cm).

L375: a suggestion – consider writing eq. (3) directly as the difference between MAR and no-MAR simulations (it is like plugging eq. (3) into eq. (4) and then Q0 and S0 are eliminated). Methodology-wise, it seems a more intuitive way to present the net recharge by MAR (and then you don’t need eq. 4).

L438-400: Please explain how you got such high infiltration rates for a field that is rated as “poor” by the SAGBI index. These infiltration rates (even the 0.06 m/hr) will redefine the site’s SAGBI rate as “good” or “excellent”. Maybe the range of Ks that you used (#1 in Table 1) was too high? You disregarded half of the failed simulations because permeability was too low for the specified upper boundary flux (L432-434). It could be that the specific site cannot accommodate the amount of applied water, and therefore you preferred the converged simulations with high Ks. In other words, it could be that the Monte Carlo approach produced an overestimation of Ks and infiltration rates?

L453: I suggest to rephrase “..accommodate more infiltration..”  (maybe to “these cells have higher infiltration capacity..”)

L462-464: Interesting result. I suggest emphasizing that your spring flooding is probably a riskier approach, as most growers will prefer to apply winter flooding when almonds (or other perennial crops) are dormant.

L469: Fig. 6 and S8 are very nice. If doable, I suggest incorporating them both as one figure in the main text.

L497-502: The dichotomic division to fine and coarse facies is arbitrary and subjective. Looking at the continuous parameters of the hydraulic functions, it is clear that all outcomes (obj. 1-3) are most sensitive to hydraulic conductivity, porosity, and residual saturation. So why divide the global sensitivity analysis into different facies?

L566-568: I would be cautious with declaring that the 3 outcomes are not sensitive to flooding frequency (especially root zone residence time). I think the source of this conclusion is site-specific – in your case infiltration rates are relatively high compared to the applied water volume. Inundation frequency of several wetting-drying cycles is common practice in MAR operations to maintain high infiltration rates (as you also stated in section 4.2.1). In SAT operations (and probably also in Ag-MAR) these cycles have also an important role in soil re-aeration.

Technical corrections

L59: change Harrington et al. (2014) to (Harrington et al., 2014)

L139: Consider changing the precipitation units to mm (60 mm and 450 mm)

L172: eq (1) – consider changing z to other letters, because you perform 1D vertical interpolation/extrapolation and in the following paragraph z denotes both resistivity [z(u)] and direction [z], which can be confusing.

L319: consider changing “domain” with simulations or realizations (as you have only one model domain).

L426: Fig. S5 – while found in the SM, the axis title in Fig. S5c should be corrected to z(m); and the x=200 in the caption should be corrected to y=200.

L458: Fig. S7 in SM, correct the legend to longer than 48 h (red/orange) and less than 48 h (blue). Also correct caption – “(green)” to “(red)”

Throughout the text – I suggest replacing “pressure” with “head” or “hydraulic head” as you also present it in units of length (e.g., Fig. S8).

L493-495: either use average or median for both facies (currently average use for fine, and median for coarse).

Fig. 10: isn’t obj. 1 title “Preferential infil.” should be changed to “infiltration rate”?