Dear authors,

Thank you for writing and submitting this interesting manuscript documenting and testing the implementation of the Richards equation in LPJ-GUESS and comparing it to the bucket water hydrology used in the model before. My background is applied mathematics, so my review will naturally focus on the method section of the paper.



General comments

Overall, the description of the method and the results is clear, the method is mathematically solid and well justified, and the detailed analysis of the results provides interesting insights into soil hydrology modelling in dynamic vegetation models. Although the employed method is not new—it is taken from Ireson et al. (2023) with two slight adaptations (unequal layer thicknesses and a sink term)—the detailed analysis of the effects of using the Richards equation instead of a bucket water hydrology makes this manuscript a valuable resource for modelers in the community, thinking about implementing the Richards equation.

The method section is in general well written and mathematically clear. The description of the baseline model seems to be split in subsections 2.3.4 and subsection 2.2.2. Potentially, the generally well-organized method section could further benefit from moving all description of the baseline model (before the updates) to section 2.3.4. There are some equations and mathematical sentences that are likely erroneous, and should be corrected before publication. However, the errors can be easily corrected and the numerics set forth is still correct and uninfluenced by the errors.

The model i.e. the partial differential equation (PDE) written in equation (5), is mathematically equivalent to the standard Richards equation. It is formulated in terms of water potential and slightly changed compared to the classical formulation using the product rule on the time derivative of the soil water content to split this term into the derivative of water content w.r.t to water potential C(Psi) and the time derivative of water potential. The only issue I find with the PDE is that the location of the sink term S within the equation is likely wrong (see line specific comment). The equation is then discretized w.r.t to space, following Ireson et al. (2023). The approach of the reference is slightly generalized by considering layers with unequal thickness, which is mathematically well justified. The result is then a system of ordinary differential equations, written down in equation (6) and (7), for which different, physically meaningful, boundary conditions are applied. For the numerical integration an off-the-shelf numerical solver, that fits the application, is used.

Two simulations have been carried out, one at a fluxtower site at Dahra and another one for the whole Sudan-Sahel region with both, new and old, model versions. The results are comprehensively analyzed and the versions are compared to each other. In addition, both model results are confronted with reference data and their performance is evaluated. For the regional simulation, the model performance is clearly improved by using the Richards equation: For the surface soil moisture, the bias of the old version is almost gone in the new version and the correlation to reference data is significantly enhanced. Additionally, the transpiration bias is strongly reduced. For the Dahra site the use of the Richards equation does not improve performance, which could well be explained by site specific confounding factors, like missing livestock grazing, as you (the authors) suggest.

A very valuable insight for modelers of the community is that you demonstrate, how the sensitivity of the model outputs w.r.t to soil texture is increased, and more boundary conditions are possible with the new soil hydrology scheme. This shows how the Richards equation is more flexible and versatile w.r.t modelling different site conditions, compared to the bucket water scheme.

In the discussion, I suggest to additionally highlight as a clear improvement of the new model version, that the Richards equation results are more physically realistic, since spatial discontinuities in the solution are avoided.

Specific comments

1. Line 145 (equation (2)): It is unclear to me whether this represents total percolation that happens between the layers per day or the percolation rate. I think to make this concrete could be beneficial.
2. Line 151: I understand that this is described in more detail in section 2.3.4. Without this further explanation I wasn't able to understand this sentence. My suggestion would be to provide the details of baseflow runoff and lateral flow runoff in the baseline model here and then only describe the changes made to these flow (if there are any) in section 2.3.4.
3. Line 167 (equation (4)): The following question is likely only because of my lack of understanding of soil hydrological modelling. In the reference Campbell (1974) the formula K = K_s(θ/θ_s)^2b+2 is used. The exponent 2b+3 is only used when an interaction term is added. However, I don't see an interaction term being used here. The question is: What was the reason for choosing to use the exponent 2b+3 instead of 2b+2?
4. Line 180 (equation (5)): The sink term S appears to be misplaced in this equation. Since the sink directly adds or removes water, it should be placed outside of the brackets of the spatial derivative, while still being multiplied with 1/C(ψ). It is not the spatial derivative of the sink that contributes to the volumetric water content rate of change dθ/dt, but the sink term S directly. Thus, i suggest to correct equation (5) as:
   
   dψ/dt =1/𝐶(ψ) (d/dz (𝐾(ψ) (dψ/dz − 1) ) − 𝑆)       (5)
   
   In the spatially discretized equation (6) the sink term Et_i/Δz_i is placed correctly, so the model code based on this should function correctly.
5. Line 184: The use of the term ordinary differential equation (ODE) is incorrect, I believe. Equation (5) is a partial differential equation (PDE), since it involves a spatial derivative w.r.t the variable x in addition to the temporal derivative. An ODE will only be achieved after the spatial discretization by the method of lines. That is equation (6) is then indeed an ODE. I suggest to use the term PDE for equation (5) and mention the "method of lines" to get from equation (5) to equation (6).
6. Line 192 ("Et_i is the implementation… "): When this equation is multiplied with C(ψ_i) both sites of the equation represent the rate of change of soil water content theta for each layer. Then the term “Et_i/Δz_i” is precisely the sink term S in the corrected version of equation (5), but not “Et_i” itself. That is because equation (5) (when multiplied with C(ψ)) also has the rate of change (time derivative) of volumetric soil water content on both sides, but not the rate of change of absolute soil water content. So, S is the additionally rate of change of volumetric soil water content due to the sink.
7. Line 240: Why is the R_drain flux calculated at the end of the day and not handled analogously to the other fluxes Es, Et and Win, that is passed to the sub-daily RE integrator?
8. Line 262 ("To ensure water mass balance..."): As I understand the technique of Ireson (2023) the use of the additional equations dQ₁/dt = q₁ and in the ODE system dQ_N/dt = q_N is not to ensure water mass balance closure, but to calculate the cumulative boundary fluxes, and thus to evaluate the deviation from water mass balance closure, i.e. to evaluate equation (13).
9. Line 266 (equation (14)): To make this equation more clear, I suggest to either use an integral over the timespan [0,t_now] or to use a sum from i=0 to i=K, over q_j(t_i) where K is the number of timesteps, i.e. sum over discrete times t₁,…,t_K. The way this equation is written now, the sum seems to be taken over a continuous interval (uncountable many terms) which is at least not a standard mathematical thing to do.
10. Line 463 ("due to a large spatial variability..."): Why does the large spatial variability in the correlation lead to a lower average correlation?
11. Line 634 ("updating the soil water dynamics in LPJ-GUESS resulted ..."): Based on my interpretation of Figure 4, Figure 5, Table 1, Figure S5 for the Dahra site simulation and Figure 7 for the regional simulation compared to GLEAM, there is indeed an improved performance for the regional simulation when using the Richards equation, but not for the Dahra site simulation. An additional advantage I see for using RE is that, for both simulations (regional and Dahra) the soil moisture profile "looks" much better using the RE, that is sharp gradients without physical meaning (which I would call artefacts of the numerical method) are avoided (see e.g. Figure S7).
12. Line 643 ("Including these boundary conditions..."): Maybe add: "may therefore further improve the overall model performance [for the regional simulation]"? For the Dahra site this seems not to be the case, as I understand the explanation of line 621-623.
13. Line 696 ("This also resulted..."): I would be interested in understanding this matrix-based approach better. Is it an implicit method (like backward Euler or trapezoidal method)? Using an implicit method (based on matrices) could indeed increase computational efficiency since no restriction of the length of the timestep is needed for stability and one may be able to avoid sub-daily time stepping all together. However, using an implicit method would also mean that a nonlinear equation system would needed to be solved each timestep.
14. Line 718 ("Indeed, from our sensitivity..."):  I was unable to locate this sensitivity analysis in the paper. Does this sentence refer to Figure 14?

Technical corrections

1. Line 148 (“see Sup. Mat. S1.1”): A typographic comment: Could you write this out as in line 131 for consistency?
2. Line 234 ("The removal of water... "): This sentence could be read to imply that the water infiltration Win is not included in the ODE solver runtime. However, paragraph 2.3.5 clearly states that Win is passed to the ODE solver analogously to Et and Es. My suggestion is to add "Win" to the list of things implemented inside the ODE solver routine, if this is the case.
3. Line 525 (Figure 10): For easier readability and more consistency it would be beneficial to use dotted lines for default values, as in Figure 2. (Or change Figure 2 respectively.)
4. Line 528 ("meteoroloigical"): This is a typo.
5. Line 618 ("the the"): This is a typo.
6. Line 653 ("down to -1.4 m"): Is the minus sign in front of 1.4m a typo?

These are all my comments. Feel free to contact me if needed.

Best Regards,

David Hötten