This is a review of "Development of Inter-Grid Cell Lateral Unsaturated and Saturated Flow Model in the E3SM Land Model (v2.0)" by Qiu et al.

The authors describe the development of a saturated lateral flow parameterization for between-gridcell water movement.  The modified model is compared to a fully 3d subsurface flow model.

Derivation of equations

In general, I thought that the derivation of the moisture movement equations could be improved.  Around line 105, the authors reference Oleson et al. [2013], which describes the conservation of mass in one dimension that is used in the CLM model.  But their equation 7 is not what is found in Oleson et al., and instead is written in the general 3d form.  The authors here use the index n to denote different control volumes, before switching to an index k used to describe the layers of the 1d soil structure.  I think it would be more clear to either 1) describe the 3d equations first, then show how the specific 1d case leads to the ELM/CLM equations described in Oleson et al., or follow Oleson et al. and then show how the 3d equations are used to define the new term in equation 14.

Regarding the equations describing the lateral flux, e.g. 14 and following, the areas used to convert between fluxes and volumes should be clarified to show whether they are actual surface areas, or projected areas.  The appropriate area is the area that is normal to the direction of the flux.  This is why I found the description of equation 16 confusing.  I don't think describing the fluxes relative to z' (the rotated z axis) make sense.  This is still a 1d column model, with the nodes aligned verticaly, therefore the fluxes in the column are all in the vertical.  Presumably the coupling to the atmosphere also occurs  in the vertical.  It would be more clear to me to note that a cosine arises in equation 16 due to the projected area of the surface being smaller than the surface area by a cosine factor.

I would also like to know if the presence of z in equation 15 is correct.  I understood the gravity term in the modified Darcy equation to be the sin(theta) term.

Why does equation 17 not have a similar form as equation 15?  Also, I would like to see the calculation of the transmissivity T described here rather than simply referenced.

What is the size of the contribution of the unsaturated lateral flow term, for both the benchmarking and application simulations?  It would be useful to indicate the relative importance of this term, and whether this impacts the simulations significantly.  It would seem to be straightforward to turn off this flux to test this issue.

Evaluation over lww

A comparison to a LWW PFLOTRAN simulation would have been interesting.  Given that PFLOTRAN was used in the benchmarking section, why was it not used in the evaluation section?

The WTD map in figure 8 shows that the addition of lateral flow helps to better resolve the uphill/downhill differences in water table depth, but there are still significant differences relative to the Fan WTD map.   The authors note that calibration of the f_d parameter to give a better match to the Fan WTD map may not be fruitful due to differences in climate forcing.  But given the relatively large differences, it would be informative to do a sensitivity test for the f_d parameter.  For example, is there a value of f_d that further lowers the water table depth, and better resolves the riparian areas as shown in the Fan map? 

Statements such as "The effects of WTD changes on the energy fluxes were more pronounced at low elevation cells, especially at the stream and its surrounding cells. The delivery of the groundwater through the lateral flow to the valleys supported higher LH while reducing the SH compared with ELMv2.0 which has little spatial WTD variations" do not appear to be well supported by figure 9.  Instead of highlighting the differences between the uphill and downhill areas apparent in figure 8, figure 9 shows spatial patterns having broad domain-wide patterns.   Why do the water table patterns in figure 8 show much more structure?  For example, larger LH values do not appear consistently in the riparian areas.  Similarly, the patterns in the difference maps only show scattered points rather than a clear riparian pattern.  Why is this?  

The comparisons to observations (figures 10 and 11) do not seem to add much insight into the relative model behaviors.  Given the authors choice to not calibrate the models, I don't think it can be stated that the differences between the observations are due to model structure.  For example, the A121 differences in figure 10 might be smaller if the f_d parameter in ELMv2.0 model had been calibrated.  The statement that both models "were able to capture the major fluctuations and wetting/drying cycles of soil moisture (SM)" seems over-stated.  The rain events are generally captured, but the magnitude of the reponse, and the dry-down rate is generally poor.  What information are the authors trying to give to the reader with these figures?  Similarly for figure 12; one does not need to perform a model simulation to be aware that shallower water tables will typically have colder temperatures, higher LH, and lower SH than deeper water tables.  Any two model versions having different water table depths would presumably show this behavior.  I don't see that this figure adds any additional insight to the results.  

The authors have the original idea of considering horizontal groundwater flow in the source/sink terms of the equation. While these are interesting results, I believe the paper could be improved by considering the following points.

About the equation

To clarify the difference between ELMlat and PFLOTRAN, it should be noted that the unsaturated horizontal flow is the results of one previous step.

Equation 16 makes sense if flux includes the surface area of the terrain, but I can't determine that with the current description.

Equation 17 is correct in the equation itself, but does it not have to account for slope gradient as in equation 15? Need a description of why unsaturated horizontal flow is considered but saturation is not.

Regarding the model benchmark

The authors mentioned that the moisture retention curves used in ELM and PFLOTRAN are different, but why not show how much they differ? Showing the results of the fitting would be helpful for the discussion.

Regarding the results

What is the reason for using top 5 layers in the comparison with PFLOTRAN as far as Figure 5 is concerned, I think it would be a fair assessment to include up to about top 10 layers. Also, what is the reason for using MAE? It tend to be small because the denominator is larger than the numerator. Why not present other indicators as well?

I understand that Figures 6 and 7 are ELMlat-PFLOTRAN. If so, some discussion of the CH and DH soil moisture and groundwater table results is needed; ELMlat has less soil moisture on the down hill and more on the up hill than PFLOTRAN, even though the groundwater table is too flat compared to PFLOTRAN. Is this solely due to the water retention curve? Need to describe the difference in equations and the effect of the solution method.

Regarding Figure 9, from Figure 8, the groundwater table is deeper near the watershed boundary, but even in such locations, does the soil temperature decrease, LH increase, and SH decrease? Are these results consistent? Results and discussion on whether horizontal unsaturated flow has an effect and how the seasonal planar distribution varies are needed.

Minor comments

Line 127 may be clearer in Figure 1 than in Figure 2

I think (d) in Figure 10 is a mistake for A149.

The equation numbers are not bracketed: lines 106, 116, 126, 128, 143.

FigureA1 should be the result of all layers, not top 5 layers.

FigureA2 does not have results for (d),(e),(f).

FigureA3 has all figures marked as (a).

Line 540 is Soil MoistureâBased (I noticed this by accident)

This manuscript gives a good description of the parameterization scheme of Inter-Grid Cell Lateral Unsaturated and Saturated Flow Model, and evaluates the simulation performance of this scheme. But I think there are still some things that need to be improved to make it better. My detailed comments are as follows:

About the parametric scheme

1. Line 94. For soil water movement, k is an important parameter. Is the description of equation 3a accurate? I noticed that ELM originated from CESM1.3, and the calculation of k in CESM1.3 is divided into two parts according to soil layer (for details, please refer to Oleson et al., 2013, equation 7.89). Why did this paper modify this? Will such modification have a significant impact on the results?
2. Line 100. The modification of equation 4 can be confusing. I recommend the author can cite “Yan Yu, Zhenghui Xie, XubinZeng, Impacts of modified Richards equation on climate simulation in the regional climate model RegCM4, Geophys. Res. Atmos., 119, 12,642–12,659, doi:10.1002/2014JD021872, 2014.”. It would be more appropriate to explain the reasons for the changes.
3. Line 130. Equation 14, does the direct modification of the tridiagonal equation mean that soil water lateral flow is carried out firstly and then vertical flow is carried out in 3D soil water flow? In fact, however, they happen at the same time. Will this have a significant impact on the results?
4. Line 137. I am confused about the expression of equation 15. Generally speaking, the soil water lateral flow is driven by soil matrix potential, but the gravitational potential (z) is obviously included in equation 15. Clarify please.
5. Line 137. I think the parametric scheme of soil water lateral flow is the focus of this paper, but this part is too simple in this paper. It is important to give a detailed process of the derivation of equations rather than simply referenced. Moreover, it seems that the final parametric equation is not given in equation 14. I suggest: 1) First introduce the soil flow scheme in ELM; 2) Then introduce the soil water lateral flow scheme and your improvement; 3) Finally, the tridiagonal equation including soil water lateral flow is given. Or switch step 1 and step 2.
6. Line 149. The description of equation 17 is simplistic. In fact, in my opinion, equation 17 is just as important as equation 15, so the derivation of 17 should be detailed rather than simply quoted, because I have noticed there are differences between equation 17 and the equation in Fan (2007).

About model benchmarking

1. In the benchmarking model, I am more interested in the comparison among ELM, ELMlat and PFLOTRAN. This comparison may be important to reveal the importance of soil water lateral flow. and can more directly show the improvement of this paper.
2. Is the lateral flow of unsaturated and saturated soil water equally important? Or is one of them dominant? How much do they contribute to soil water movement?

About LWW

1. Line 297. Why did this paper focus on the summer soil temperature?
2. For LWW case, authors only verified soil moisture for shallow soils (25cm), but I think the results below 1m are just as important. I think more soil moisture data (such as ERA5-land or satellite inversion data) can be used to verify the simulated soil moisture, soil temperature and other results, rather than just using two stations data. This is important to illustrate the importance of soil water lateral flow.
3. Figure 8. In addition to the differences between Figure b and Figure c, the resolution of Figure b seems to be lower. Is this related to the modification of the depth of the soil layer? If yes, can more reasonable soil depth distribution be considered in future work?

Minor comments

1. Is line 138 and line142 duplicated?
2. Line 146. The unit of is m3/s, but the unit of q in line 124 is mm/s. There is an obvious difference between the two. Whether this is a clerical error, if not, please explain.
3. It will be better to label (a) (b) (c) in Figure 12.