Comments on: Forward and inverse modelling of water flow in unsaturated soils with discontinuous hydraulic conductivities using physics-informed neural networks with domain decomposition.

Summary 

This paper presents the results from a comprehensive study using PINNs as a forward and inverse numerical solution for the Richardson-Richards equation. They tested new approaches for applying the PINN method, including a layer-wise locally adaptative function intended to work with layered heterogeneous soil profiles. In addition, the authors compared their approach to well-known numerical solutions for the Richardson-Richards equation, namely Finite Difference and Finite Element Methods (FDM and FEM). The PINNs approach was also validated with soil moisture measurements performed in a soil column in controlled conditions.  

The paper appears to be relatively novel, being the first application of PINNs to the Richardson-Richards equation (to my knowledge). The literature proposed and the figures presented are of high quality. I enjoyed reading the paper. The results are encouraging on the applications of PINNs to model hydrodynamics in porous media, even if it takes much more time when compared with the classical approaches. We don't need to impose well-known boundaries and initial conditions, which is attractive once they are difficult to obtain in field applications. The domain decomposition for the layered soils is also very promising. Even classical approaches such as FDM and FEM struggle with heterogeneous soil profiles. So, I think the PINNs with the domain decomposition did quite well in modelling the soil water dynamics in the soil column. 

Specific comments/questions (that should be addressed and commented before publication): 

• It would be interesting to test the inverse solution with soil matric potential measurements (data is available if needed).  
• What is your opinion on going to 2 and 3D modelling? Could the domain decomposition proposed in the paper be applied to speed up 2 and 3D solutions? I think that would be the actual gain in this methodology. FEM applications for the fully 3D solution of Richardson-Richard's equation are still slow and have many complications with mesh, especially for large domains. This also applies to the boundary and initial conditions imposition. 
• What about non-Darcian conditions, macropore flow, very high clay content soils. Do you think the method could be applied?
• What about root-water-uptake? How can this be included in your approach? There exist some analytical solutions for these problems (Yuan and Lu, 2005[1])
• Do you think one day the PINNs could take over the classical approaches? What is limiting it?
• What about practical applications? Irrigation management or contaminant transport in the vadose zone.

Overall, the paper is well written. The sections are balanced, and the flow is good, making the paper enjoyable to read.

  [1] Yuan, Fasong, and Zhiming Lu. "Analytical solutions for vertical flow in unsaturated, rooted soils with variable surface fluxes." Vadose Zone Journal 4.4 (2005): 1210-1218.

Review of “Forward and inverse modeling of water flow in unsaturated soils with discontinuous hydraulic conductivities using physics-informed neural networks with domain decomposition” by Bandai and Ghezzehei.

In this manuscript the authors tested a physics-informed neural networks (PINNs) method to solve the Richardson-Richards equation for simulating unsaturated soil water dynamics. The authors also investigated the capability of the method for obtaining inverse solutions. As coupling data-driven and physics-based approaches have received much attention these days, the topic fits well with the scope of HESS. The authors have done a great job on demonstrating how PINNs performed when simulating unsaturated water flow in soils and showing applicability and limits of the method. Although the paper was well organized and written, I believe the paper has a room for some improvement. I have some comments that should be addressed prior to accepting this paper for publication. For my curiosity, I am wondering if this approach can be applied to simulate preferential type flow in soils. Is it going to be straightforward? Does it require some modifications in the model? If it can be applied to such phenomena, it would be a great breakthrough in the field of soil physics and hydrology.

General comments:

In Fig. 5, the evolution of PINNs solution is plotted. At the initialization, some of the solutions are beyond the limit of the water content as the water content values are greater than the saturated water content. Would it be possible to put some constraints to the solutions in PINNs? If so, would that improve training and overall performance? Any discussions on this matter will helpful for those who are interested in using this method. A similar question goes to the inverse solutions. I am wondering if any constraints can be applied to the target parameters that are inversely estimated. There is always a need to put some constrains to the parameters being estimated.

In the demonstration of getting inverse solution with PINNs, the authors used a 2-layered soil system. Why? If the boundary fluxes are being estimated, wouldn’t be better to start with a homogenous case? Was there a specific reason that the layered soil system was used in this demonstration?

Specific comments:

L189: It sounds a bit strange to say that soil dynamics is “controlled by the volumetric water content at the bottom.”

L193 (Eq.15): A little bit more explanations will be helpful to understand this transformation. I have no idea why the beta value gives better initial guess.

L311: If the logarithmic transformation of water potential is used, the approach is limited to “unsaturated” systems. But there are many cases you will have both positive and negative potential values. How do you deal with that?

Figure 3(b): There are some systematic differences between FDM and PINNs. Why? Are these because of the choice of spatial and temporal discretization in FDM?

Figure 10: Looks like something is wrong with the texts at the top of the figures.

Figure 16: For all three cases, the PINN solutions show that the inversely estimated initial surface flux is much smaller than the true flux. Are there any specific reasons for this?

The paper is very interesting and introduces a physics-informed neural networks (PINNs) method in a Richards’ equation context, aimed at approximating the solution to the RRE using neural networks while concurrently matching available soil moisture data. In particular, in this paper authors consider domain decomposition for handling infiltration into layered soils.

The topic is definitely up to date, the paper is well written, and provides all the details for implementing and understanding this approach. Nevertheless, I think authors should address some comments and issues before it can be accepted.

• This PINN approach appears really fascinating because it allows to integrate physics-based models (such as RRE) with machine learning features. Authors ascribe the uncertainties in Richards’equation to the choice of boundary conditions, which is surely right. Nevertheless, I think they do not consider the (even more) cumbersome uncertainties arising from the choice of model parameters, which are the result of some non-linear fitting in laboratory experiments (I am referring to the parameters in the WRC and HCF). This point is a main concern for me: as a matter of facts, unsaturated flow dynamics strongly relies on functions parameters, rather than on ICs and BCs, which are generally easier to assess. On the other hand, I see authors already published a paper on this topic: I think it would be valuable to stress the differences between the two papers
• Lines 197-198: few more words for sketching how the partial derivatives are computed would be valuable
• Figure 1: I think there is a typo in the box “Physics and Data Constraints”, since the partial derivative at the left-hand side should be accomplished with respect to time.
• I understand that the residual is computed between the synthetic data and the computed (by the PINNs) ones; in this framework, what is the rationale of comparing the PINNs output with any Richards solver (as Hydrus)?
• As far as I understand, the power of this approach is to combine physics-based models with data driven ones; according to my knowledge, this is also the spirit of Data Assimilation (DA) methods, which incorporate measurements into a physics based model, albeit  in a very different framework; these methods have also been treated in Richards’ equation context (see for instance Berardi et al CPC https://doi.org/10.1016/j.cpc.2016.07.025, Medina et al HESS https://doi.org/10.5194/hess-18-2521-2014, Liu et al JoH https://doi.org/10.1016/j.jhydrol.2020.125210 ); what is authors’ opinion about this? What would be the pros and cons of PINNs approach with respect to DA one?  Also DA methods allow to assimilate boundary conditions, as in this case, and hydraulic parameters, as well as states.  As a matter of fact, with respect to DA methods, this PINNs approach seems to me more on the theoretical side (which is definitely fine, of course) rather than application oriented.

• Authors mention the possibility to drop loss terms for IC or BC at line 225. However, they have not presented any experiment for this scenario. Could you please comment on this ill-posed configuration? How would it perform with respect to classical solver?
• Figure 11 and 3. Please replace “Fintie” with “Finite”.
• Authors make use of synthetic data: I had hard times to find where the reference to used data is described. Of course the use of synthetic data is fine, but they should highlight it at the beginning of the paper. Moreover, could you please explain how your method of synthetic data generation could compare to real measurement data? In other words, how robust is your result with respect to outliers, sensor noise and other technical issues when it comes to real data?