It rains and then? Numerical challenges with the 1D Richards equation in kilometerresolution land surface modelling
 ^{1}Institute for Atmospheric and Climate Science, ETH Zurich, Zurich, Switzerland
 ^{2}Deutscher Wetterdienst DWD, Offenbach, Germany
 ^{a}now at: Federal Ofﬁce for Meteorology and Climatology, MeteoSwiss, Zurich Airport, Switzerland
 ^{1}Institute for Atmospheric and Climate Science, ETH Zurich, Zurich, Switzerland
 ^{2}Deutscher Wetterdienst DWD, Offenbach, Germany
 ^{a}now at: Federal Ofﬁce for Meteorology and Climatology, MeteoSwiss, Zurich Airport, Switzerland
Abstract. Over the last decade kilometerscale weather predictions and climate projections have become established. Thereby both the representation of atmospheric processes, as well as landsurface processes need adaptions to the higherresolution. Soil moisture is a critical variable for determining the exchange of water and energy between the atmosphere and the land surface on hourly to seasonal time scales, and a poor representation of soil processes will eventually feed back on the simulation quality of the atmosphere. Especially the partitioning between infiltration and surface runoff will feed back on the hydrological cycle. Several aspects of the coupled system are affected by a shift to kilometerscale, convectionpermitting models. First of all, the precipitationintensity distribution changes to more intense events. Second, the timestep of the numerical integration becomes smaller. The aim of this study is to investigate the numerical convergence of the onedimensional Richards Equation with respect to the soil hydraulic model, vertical layer thickness and timestep during the infiltration process. Both regular and nonregular (unequally spaced) grids typical in land surface modelling are considered, using a conventional semiimplicit vertical discretization. For regular grids, results from a highly idealized experiment on the infiltration process show poor numerical convergence for layer thicknesses larger than approximately 5 cm and for time steps greater than 40 s, irrespective of the soil hydraulic model. The velocity of the wetting front decreases systematically with increasing time step and decreasing vertical resolution. For nonregular grids, a new discretization based on a coordinate transform is introduced. In contrast to simpler vertical discretizations, it is able to represent the solution secondorder accurate. The results for nonregular grids are qualitatively similar, as a fast increase in layer thickness with depth is equivalent to a lower vertical resolution. It is argued that the sharp gradients in soil moisture around the propagating wetting front must be resolved properly in order to achieve an acceptable numerical convergence of the Richards Equation. Furthermore, it is shown that the observed poor numerical convergence translates directly into a poor convergence of infiltrationrunoff partitioning for precipitation time series characteristic of weather and climate models. As a consequence, soil simulations with low resolution in space and time may produce almost twice the amount of surface runoff within 24 hours than their highresolution counterparts. Our analysis indicates that the problem is particularly pronounced for kilometerresolution models.
Daniel Regenass et al.
Status: closed

RC1: 'Comment on hess2021426', Jan Hopmans, 30 Sep 2021
Initially, when accepting to review this manuscipt, I was very interested to learn about more recent advances in numerical modeling. However, it took me some time to realize that the manuscript assumes knowledge at the landscape scale, and is seeking to describe LS modeling to the km scale, at much higher spatial and temporal resolution. This reviewer's expertise is in soil and vadose zone hydrology which simulates subsurface water flow using the Richards equation at the cm (m) and hour (minute) scales.
I do fully understand the inherent compiexities when taking the landscape scale as the benchmark and using those experiences to address unsaturated water flow at much higher spatial and temporal resolutions. However, over the past few decades, the soil/vadose zone hydrology field has largely advanced, and has developed computer simulation algorithms that can address these high resolutions with satisfactory computer times including multidimensional flow at increasing space (time) scales (more than cm/minute scales).
It was therefore surprising to learn that the authors have largely neglected to reference that work, and ignored those more recent developments of the past 1020 years in the soils and vadose zone hydrology literature. Instead, this manuscript makes simplistic assumptions that are no longer necessary, such as considering multidimensionality through an empirical sink term, and simple boundary conditions such as zero flux lower boundary and a nondynamic infiltration term.
Just for the author’s information, I list hereby just a few references that they could read (from a much larger list of additional references), and these are:
Recent Developments and Applications of the HYDRUS Computer Software Packages
JiÅí ŠimÅ¯nek,Martinus Th. van Genuchten,Miroslav Šejna,at https://acsess.onlinelibrary.wiley.com/doi/full/10.2136/vzj2016.04.0033
Adaptive Grid Refinement in Numerical Models for Water Flow and Chemical Transport in Soil: A Review R. S. Mansell,Liwang Ma,L. R. Ahuja,S. A. Bloom,
https://doi.org/10.2136/vzj2002.2220 at: https://acsess.onlinelibrary.wiley.com/doi/abs/10.2136/vzj2002.2220
Sustainability of irrigated agriculture in the San Joaquin Valley, California
Gerrit Schoups, Jan W. Hopmans, Chuck A. Young, Jasper A. Vrugt, Wesley W. Wallender, Ken K. Tanji, and Sorab Panday at: https://www.pnas.org/content/102/43/15352
And
Modeling Soil Processes: Review, Key Challenges, and New Perspectives
Vereecken,* A. Schnepf, J.W. Hopmans, M. Javaux, D. Or, T. Roose,et al. Modeling Soil Processes: Review, Key Challenges, and New PerspectivesH. Vereecken,* A. Schnepf, J.W. Hopmans, M. Javaux, D. Or, T. Roose et al. 2016. at: https://acsess.onlinelibrary.wiley.com/doi/epdf/10.2136/vzj2015.09.0131

AC1: 'Reply on RC1', Daniel Regenass, 30 Nov 2021
We thank Jan Hopmans for his review and the suggested literature on modelling vadose zone processes. In particular, we read with great interest the article on recent developments in the HYDRUS software package by Šimůnek et al. (2016) as it gives us a good insight on what is considered stateoftheart in the vadose zone community.
We believe that the reviewer (along with reviewer 3) might have misunderstood the aim of this paper and we would greatly appreciate feedback on how to reformulate the manuscript to make its aim more clear. It is not the intent of the authors to suggest a novel numerical method to solve the Richards equation. In particular, the solver presented here has already been introduced by Schlemmer et al. (2018) and it is based on the solver used in the TERRA ML land surface model (LSM) (Heise, 2003) and we hope that this has been made transparent in sections 1 and 2 of our paper. The LSM community and (as part of it) the authors are aware of the progress that has been made in the numerical solution of the Richards equation in the last twenty years. As an example, we explicitly mentioned the recent implementation presented in Lawrence et al. (2019) based on the method introduced by Kavetski et al. (2001), which is comparable to what is used in the HYDRUS model (Šimůnek, 2009). We also referenced publications that address a threedimensional implementation (Maxwell et al. 2015, Mastrotheodoros et al. 2020).
Note that the current version of the CLM model (Lawrence et al., 2019) is a noteworthy exception. In its most recent implementation, the numerical solution of the Richards equation uses an adaptive time step. In contrast, virtually all European centers and modelling consortia use a solver that is similar to the one found in the TERRA ML LSM (Heise, 2003), and use the same time step as in the atmospheric model. This paper is an attempt to raise awareness in the weather and climate modelling community that the numerical methods in the corresponding LSMs are often insufficient to achieve numerical convergence of the Richards equation given high precipitation intensities. This is an urgent matter, as higher precipitation intensities are a natural consequence of the transition in weather and climate modelling to the kilometer scale (see e.g. Ban et al., 2014).
Clearly, LSMs should be adapted to face this challenge and one evident way forward would be to 'catch up' with the progress made in vadose zone physics in the last twenty years. While this is certainly a reasonable and valid approach, we will in the following take the HYDRUS 1D model as an example to highlight potential issues that would arise when using a comparable model as a module inside a weather/ climate model. Eq. 7.1 in Šimůnek et al. (2009) is solved with a Picard iterative solution scheme. Following the notes on discretization (Šimůnek, 2021) this means in practice that up to 10 iterations per time step are required to achieve convergence in the presence of sharp gradients in hydraulic head (which is the case during the infiltration process, in particular during convective events in summer). Moreover, again following the notes on discretization (Šimůnek, 2021), in the presence of sharp gradients, the vertical discretization should be in the order of cm and the 'stretch factor' between two vertical layers should not exceed 1.5. Note that these recommendations are pretty much in line with the findings presented for the Schlemmer et al. (2018) solver investigated in our manuscript in a clear and systematic manner  with the exception that a stretch factor of 1.5 would probably be too high for the scenarios presented, as it would roughly correspond to the grid layout with 11 layers. We are concerned with the application of such a solver in fully coupled weather and climate models. For convection resolving simulations, implementing an iterative solver for the Richards equation would on average probably increase the time to solution of the hydrology module by a factor of three to five at most as the time steps in these models are usually shorter than one minute. Compared with the overall time to solution of the fully coupled model this would probably be acceptable but not optimal as variations in the time to solution are largely unwanted in timecritical applications such as numerical weather prediction, especially in operational contexts. More issues would arise due to the higher vertical resolution as this would not only increase time to solution but also the memory footprint of the model. Given our results and the fact that most soil models have soils deeper than 2 m, it seems unlikely to achieve an acceptable convergence with less than 3040 vertical layers. As highlighted in the manuscript, this is roughly a factor 34 more than what most LSMs use as a standard setup. Given that, it is uncertain, whether the Richards equation captures the relevant physics on the scales of interest (Beven and Cloke, 2012) it is not evident whether the weather/ climate modelling community would be willing to invest that amount of resources on the solution of the Richards equation. In this context, our manuscript could serve the weather/ climate modelling community as a starting point for a discussion on potential ways forward.It is also worth pointing out that the same modeling framework is often used for horizontal resolutions ranging from O(100 km) to O(1 km), increasingly also in global climate models (Stevens et al. 2020). This wide range of resolutions implies significant challenges both on conceptual and computational levels.
References
Ban, N., Schmidli, J., and Schär, C. (2014). Evaluation of the convectionresolving regional climate modeling approach in decadelong simulations. Journal of Geophysical Research: Atmospheres, 119(13):7889–7907.
Beven, K. J. and Cloke, H. L. (2012). Comment on “Hyperresolution global land surfacemodeling: Meeting a grand challenge for monitoring Earth’s terrestrial water” by Eric F. Wood et al. Water Resources Research, 48(1).
Heise, E., Lange, M., Ritter, B., and Schrodin, R. (2003). Improvement and validation of the multilayer soil model. COSMO newsletter, 3:198–203.
Kavetski, D., Binning, P., and Sloan, S. (2001). Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of Richards equation. Advances in Water Resources, 24(6):595–605.
Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W., Swenson, S. C., Bonan, G.,Collier, N., Ghimire, B., van Kampenhout, L., Kennedy, D., et al. (2019). The Community Land Model version 5: Description of new features, benchmarking, and impact of forcing uncertainty. Journal of Advances in Modeling Earth Systems, 11(12):4245–4287.
Schlemmer, L., Schär, C., Lüthi, D., & Strebel, L. (2018). A groundwater and runoff formulation for weather and climate models. Journal of Advances in Modeling Earth Systems, 10(8), 18091832.
Šimůnek, J., Van Genuchten, M. T., & Šejna, M. (2016). Recent developments and applications of the HYDRUS computer software packages. Vadose Zone Journal, 15(7), vzj201604.
Šimůnek, J., Van Genuchten, M. T., & Sejna, M. (2009). The HYDRUS1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variablysaturated media. University of CaliforniaRiverside Research Reports, Version 4.08.
Šimůnek J. (unknown date). Notes on Spatial and Temporal Discretization (when working with HYDRUS). Retrieved from http://www.pcprogress.com/Documents/Notes_on_Spatial_and_Temporal_Discretization.pdf, on November 14, 2021
Stevens, B. and et al., 2020: The Added Value of Largeeddy and Stormresolving Models for Simulating Clouds and Precipitation. J. Meteorol. Soc. Japan, 98 (2) , 395435 https://doi.org/10.2151/jmsj.2020021

AC1: 'Reply on RC1', Daniel Regenass, 30 Nov 2021

RC2: 'Comment on hess2021426', Anonymous Referee #2, 08 Oct 2021
The comment was uploaded in the form of a supplement: https://hess.copernicus.org/preprints/hess2021426/hess2021426RC2supplement.pdf

AC3: 'Reply on RC2', Daniel Regenass, 09 Dec 2021
We thank the reviewer for their careful review and their comments on our manuscript. In the following, we will discuss their major comments:
 The first major comment addresses the use of the saturationform of the Richards equation (as used in our paper), as opposed to the headbased formulation. We fully agree with the reviewer that the mixedform of Richards’ equation is to be preferred over the saturation form – mainly because of the issues pointed out in the appendix. We are currently working on an implementation of the mixed form and plan to revisit our analysis once the new solver is ready. It is however important to stress that most European weather centers/ weather and climate modelling consortia use the saturation form of Richards’ equation (see manuscript lines 147150). The aim of this paper is not to suggest a novel approach to the numerical solution of the Richards equation, but rather to raise awareness on the (poor) convergence of Richards’ equation as used in weather and climate models and to shed light on the associated errors. Concerning the upper boundary condition (i.e. the formulation of infiltration): We read the corresponding formulation (Eq. 21) in Tubini and Rigon (2021) and think that our numerical implementation is adequate in the saturation form (Eq. 9) as it is a flux boundary condition prescribed by the precipitation flux. Note however, that in the saturation form, an upper limit for this flux must be prescribed as the available pore volume is finite and the flux cannot exceed saturated hydraulic conductivity, because the suction (diffusive) term must vanish at saturation.
 The term ‘successfully’ might indeed be misleading. Here we mean that land surface models using the saturation form of Richards’ equation are used in land surface models and yield physically meaningful results. This does however not imply that the mathematical/ numerical implementation is correct in a rigorous sense. We will try to clarify this issue in the revised version of the manuscript.
References
Tubini, N., & Rigon, R. (2021). Implementing the Water, Heat and Transport model in GEOframe WHETGEO1D v. 1.0: algorithms, informatics, design patterns, open science features, and 1D deployment. Geoscientific Model Development Discussions, 147

AC3: 'Reply on RC2', Daniel Regenass, 09 Dec 2021

RC3: 'Comment on hess2021426', Anonymous Referee #3, 21 Oct 2021
The authors presented in the manuscript a study dealing on ‘improvements’ in temporal and (vertical) spatial resolution in LSM. From the title alone I was attracted as the title suggest novelty and analysis, whereas the manuscript itself does not show any novelty at all. For me it seems y that some basic soil physical principles are not totally clear to the authors. Additionally, recently, Vereecken et al. (2019) published an overview article (Infiltration from the Pedon to Global Grid Scales: An Overview and Outlook for Land Surface Modeling), where from page 8 onwards exactly the problems raised here are shortly discussed. Additionally, there is a link provided (https://www.pcprogress.com/Downloads/Public_Lib_H1D/Using_HYDRUS1D_to_Simulate_Infiltration.pdf) to a full set of simulations using a physical based model showing the effects of different spatial and temporal discretization on infiltration.
Additionally, there are many flaws and sometimes easy problems are made complicated. I worked over the manuscript in detail to about half and made many suggestions and provide comments either in the attached scan or below. After, I decided to have only a quick look at the results as I already made my opinion that the manuscript in its current form is not appropriate for publication.
ABSTRACT
Line 5: If you want to represent runoff and infiltration correctly, even subhourly data are needed, as short but intensive rainfall (e.g., thunderstorms) often cause high runoff as they often reach dry soils.
Last sentence in the abstract: I do not agree, as this is fact irrespectively of the scale models run at (as I also believe that the models commonly used will not change in their algorithms if smaller horizontal resolution will be used). To my understanding the scale (horizontal resolution) only changes the driving forces, especially precipitation.
INTRODUCTION
Line 63: This is exactly the problem in LSM. Here, it was suggested to use a depthdependent profile of saturated hydraulic conductivities. Maybe by this more water can infiltrate but still we will have a wrong model setup generating good results. For sure, hydraulic conductivity will change over the depth profile (which has been shown in many publications and often Ks is high in ploughed layers) but why not first take care that the infiltration process is described correctly?
Line 7580: To disentangle the effect of how Ks will be defined at different scales and what the effects of different forcing are (here precipitation) a simple scaling exercise could have been performed, such as those done for crop models and water balances by Constantin et al. (2019, doi.org/10.1016/j.agrformet.2019.05.013)
Line 82: by some form of the Richards equation. Exactly, this is the problem, as the Richards equation will be solved in the diffusivity form in most LSM and not in the mixed form as most (small) scale soil hydrological models. And this exactly generated the problems as an ‘external’ infiltration concept has to be implemented to solve the infiltration/runoff partitioning (such as GreenAmpt for example). An overview that this holds for most LSM is provided in Tab. 2 of the Vereecken et al. (2019) already mentioned above.
Line 8687. Indeed, different PTF will allow higher or lower infiltration but with respect to the problem analysed here this is not a solution. Would be the same ‘wrong change’ as already done for the vertical dependent Ks.
Line 8892. Indeed, if high resolution maps are used it does not make sense to use class PTFs for e.g. the 12 (or in some LSM only 11) USDA soil classes, as the same output will be generated as we would have for coarser resolution of the soil maps.
Line 9293: Indeed, as from a soil hydrological perspective this has been already analysed over the last 2030 years but seems not to be percolated to the LSM community (or at least part of this community).
Line 100: Adaptive timestepping. At least in the soil hydrological model I use, adaptive time stepping has been implemented since the early 1990s.
Overall, it would be good just to benchmark the LSMs against analytical solutions of infiltration or against soil hydrological models and one would see if they can capture infiltrationrunoff partitioning, infiltration front propagation etc.
THE SYSTEM UNDER CONSIDERATION
This section is somehow textbook knowledge and remind me of what I show in my lectures!
Line 121123: The features you encountered in the experiment are not necessarily connected to soil heterogeneities. It is known that even homogeneous sandy soils will show so called finger flow, which are caused by the hydrophobicity of the sand grains. One a grain will be wetted, the sand turns to be hydrophilic and water can enter the profile, while other regions will stay dry. Additionally, this section is a bit out of scope here and mixes up different processes and problems with those you are working on.
SOIL HYDRAULIC MODEL
Line 159: Here, you state that the two soil hydraulic models need parameterization, which depend on soil texture. This is not fully correct. They depend on the soil pore (distribution), which depend on one hand on soil texture, but also on soil structure. For the conductivity function it is even more complicated as it also depends on pore connectivity.
Line 159. Where did you get the parameters from for the loamy soil? Which PTF or better class PTF was used? Looking at the parameters presented in Tab. 1 I wonder why you showed m for MvG and not n as m is classicaly bound to n by M = 11/n. I also wonder about the alpha value which is 3.6 and therefore too large for classical loamy soils. Finally, why is the saturated hydraulic conductivity Ks different between the two models used as Ks is a material constant of the soil itself?
GROUND RUNOFF
Maybe I got something totally wrong here but isn’t this Dunne overland flow if the soil profile is saturated (often also denoted as saturation excess overland flow)? If yes, this is totally different (in its physical principles) as runoff generated by Hortonion overland flow, which is the precipitation excess with respect to the infiltration capacity?
Finally, why assuming a zeroflux boundary condition at the lower end. To my understanding and neglecting lateral water flow in the saturated zone this will later or earlier always generate Dunne overland flow if the ET is smaller as precipitation.
BOUNDRAY CONDITION
Here we are. In your case, Imax is set to Ks. In reality, this is only the case if infiltration will be only driven by gravity and not by capillary forces. This might happen if the soil is close to saturation but if the soil is much dryer infiltration at its onset can be larger as Ks. This is simple infiltration theory.
NUMERICAL IMPLEMENTATION FOR EQUIDISTANT GRIDS
Line 196: Are you sure that most LSM solve the Richards equation on equidistant discretization. I do not believe so, as to my knowledge most have finer discretization at the top end and increase the spacing with depth. This is also stated in line 216, where it is mentioned that LSM commonly choose grids with finer spacing close to the surface.
Finally, I do not understand why you have to stress on implementing nonequidistant gridding in your model as this is stateoftheart in all other models.
Line 199: There are two points mixed up. The special and temporal resolution.
EXPERIMENT A…
Line 273: …reproduce the analytical solution. What is the analytical solution? Physically, maximum infiltration at saturated conditions without ponding water at the surface equals Ks.
To not understand what you mean by convergence. Numerical convergence does not necessarily mean that the process is represented physically (e.g. shape and wetting front speed)
EXPERIMENT B
Actually, it is not really of interest here how the forcing was calculated by COSMO as any forcing could have been used to show what you want to show. I would either delete or shift to the Annex (if any Annex is allowed in HESS)
FIGURES
Figure 5: Sorry if you define Theta_f as volumetric water content as done in Tab. 1 how can you reach Theta_f of 1 here. This would mean only water and no soil in the systems! I would denote this as lake or ocean!
TABLES
In general, headers for tables should be above the table and not below (see Tab 1 and 2).

AC2: 'Reply on RC3', Daniel Regenass, 09 Dec 2021
We thank the reviewer for taking their time to review our manuscript. First, we would like to stress that the scope of this paper is not on ‘improvements’ but rather to systematically asses numerical issues with the Richards equation as implemented in today’s stateoftheart land surface models (LSMs) coupled to weather and climate prediction models, and to highlight their interaction with increasing precipitation intensities. Since both reviewer 1 and 3 seem to have understood our manuscript as a suggestion for a new solver, we acknowledge that the actual scope has to be made clearer in the introduction. Any feedback on how to make the scope and aim of this paper more understandable would be greatly appreciated. It is clear to the authors that the stateoftheart in many if not most land surface models (e.g. Heise et al., 2003, Balsamo et al., 2008) is vastly different from the stateoftheart in infiltration models used in the vadose zone community (e.g. Šimůnek et al, 2016). Yet, it must be kept in mind that these LSMs are run under operational constraints, with severe limitations on the execution times. This is well documented in the review by Vereecken et al. (2019) and we thank the reviewer for pointing us to this review. While Vereecken et al. (2019) touch on the subject of numerical implementation, they do not discuss numerical convergence analyses, and nor the role of such analyses for establishing the magnitude of respective numerical errors. We thus still believe that systematic convergence analyses of the Richards equation yield added value for the discussion on how to treat the infiltration problem in LSMs. Moreover, the interaction of numerical errors in the solution to Richards equation with precipitation intensities in the context of weather and climate modelling has so far received little attention and are this topic is thus novel for the LSM community. Systematic analyses of numerical errors give us a better baseline for the decision, what to spend the limited execution time on.
In the following, we will address the individual points raised by the reviewer. The reviewer also mentions flaws that they pointed out in handwritten annotations in the manuscript (see supplement to RC3). There is a considerable overlap between the points raised directly in the review and the handwritten annotations in the supplement. We will thus try to address the handwritten comments, where they are not typos/ formal issues and where we feel that they do not correspond to a point raised in the review and mark the corresponding answers with (IA) for ‘inline annotation’.
Line 5: Yes but the feedback to the atmosphere will take place on hourly to seasonal time scales. This is what the weather/ climate modelling community is mostly concerned with.
Line 10 (IA): Here we mean grid spacing.
Line 17 (IA): We think that this will become clear from the experiments, but we could say fast increase in layer thickness in this case means exponential increase in layer thickness as is the case for geometric grids (Table 3).
Line 2123: The solver in the LSM will not change when changing the horizontal resolution of the LSM. But in a coupled setup, the LSM will be subject to a different precipitation forcing. This stems from (a) the effect of smoothing of the simulated precipitation in a coarser atmospheric model and (b) the effect of the parameterization of convection. While parametrized convection in largescale models tends to produce more drizzle and less highintensity precipitation, in kilometerscale models the convective precipitation is simulated explicitly leading to more highintensity precipitation. Both effects are further explained in the introduction and disentangled in chapter 4.3.
Line 4854 (IA): The equations of a dynamical core in the atmosphere are strictly based on conservation laws (mass, momentum, energy) that are described by some simplification of the NavierStokes equations. We do not have an according set of equations to describe conservation laws for the land surface. Therefore, land surface modelling (along with other subgridscale parameterizations in weather and climate models) is subject to approaches that are more heuristic. While this is reasonable given the complexity of the matter, it makes the interactions of the equations at hand (and their numerical representations) harder to understand and quantify.
Line 63: We deliberately choose an implementation of Richards equation that is used in a weather and climate model (Schlemmer et al.,2018). Depth dependent profiles of hydraulic conductivity may be found in similar forms in other LSMs (e.g. Decharme, 2006; Ducharne, 2016). We fully agree that infiltration processes are often poorly represented in weather and climate models. We are however convinced that the solution to the governing equations must be represented accurately by the solver, and to this end convergence studies such as the one presented in our manuscript are important and should be carried out complementary to the implementation of additional processes and a better process representation.
Line 7580: Again, this is mostly a problem of the forcing precipitation data in the coupled model. Here we are specifically concerned with the infiltration process in weather and climate models, where rainfall intensities are not solely subject to spatial aggregation (i.e. the size of the grid cell or LSM tile in case the rainfall is disaggregated), but also the way precipitation is represented. Mesoscale weather models and regional climate models are often operated on the kilometerscale, where convective processes need not to be parameterized in contrast to global models. There are marked differences in precipitation intensities between socalled convection permitting (or convection resolving) models and models with parameterized convection processes (Ban et al., 2014). The scaling of hydraulic properties such as hydraulic conductivity is beyond the scope of this paper, but we agree that the subgridscale heterogeneity should urgently be addressed in future research. Thank you for pointing us to the reference.
Line 82: We fully agree with the reviewer that the mixed form of Richards equation should be used and we are currently revisiting our analysis for the mixed form. However, as can be seen in Vereecken et al. (2019), Table 2, four out of seven of the listed LSMs use the diffusivity (saturation) form of Richards’ equation. Moreover, five out of seven of the listed LSMs generate surface runoff when the precipitation rate exceeds saturated hydraulic conductivity. Therefore, a convergence analysis on such an implementation is of practical relevance to the LSM community.
Lines 8687 and 8892: There are ongoing attempts to improve the representation of hydrological processes along several different pathways. None of those will be efficient on its own. We agree with the reviewer that the infiltration process is a particularly critical one, and it must first be represented accurately by the corresponding governing equations and their numerical representation. Otherwise errors in the dataset will interact with physical and numerical errors, which hampers LSM development.
Lines 9293: The lack of awareness in the LSM community – or, more broadly speaking, the weather and climate community – is precisely the issue we want to address with our manuscript. This would allow for a discussion on ways forward. We agree with the reviewer that this ‘knowledge gap’ should be closed in order to foster the collaboration amongst hydrologists, soil physicists and LSM modelers.
Line 100: Variations in timetosolution is problematic in the operational context of weather and climate modelling. This might at least partially be the reason, why parts of the LSM community are hesitant to adopt such approaches.
Line 133 and line 150 (IA): In TERRA ML (Heise et al., 2004, Schlemmer et al., 2018), E is a sink term in the root zone. We agree with the reviewer that it would be more elegant to include transpiration in the suction head (if this is what they mean). However most LSMs are simply not that far in their evolution yet.
Line 146 (IA): Q is a sink term that should emulate horizontal transport in complex terrain to some extent. We follow the implementation of Schlemmer et al. (2018). While this approach is of course very heuristic, it keeps the LSM onedimensional and it is computationally efficient.
Benchmarking LSMs against infiltration models is an excellent idea and we are very open to this suggestion.
Subsection 2.2.1 “The system under consideration”: We are open to omit this section. We believe however that it nicely illustrates the problem with vertical resolution. If one compares the sharp moisture gradient across the wetting front with the vertical discretization in stateoftheart LSMs, the issue as well as the effect on infiltration should be immediately clear.
Subsection 2.2.2 “Soil hydraulic model” (line 159): The values are adopted from the ORCHIDEE LSM (Ducharne, 2016). We are open to revisting the analysis with different values, but it will likely not change the results of the convergence analysis tremendously.
Subsection 2.2.3 “Ground Runoff”: The parameterization of ground runoff in our version of TERRA ML (Heise et al., 2003) is described in Schlemmer et al. (2018). In contrast, both saturation excess and infiltration excess are included in the surface runoff. In our case it would correspond to a saturation of the upper soil and an exceedance of saturated hydraulic conductivity respectively.
Subsection 2.2.4 “Boundary condition”: The answer to this point is partially included in the answer to the point raised concerning using the mixed form of Richards’ equation (line 82). As pointed out above, five out of seven of the listed LSMs in Table 2 of the review by Vereecken et al. (2019) generate surface runoff when the precipitation rate exceeds saturated hydraulic conductivity. While this approach is likely not valid for very dry soils, it is still of practical relevance as this assumption is included in many if not most LSMs and it is probably also a fairly good approximation for wet and moderately dry soils (see also Ogden et al., 2017). However, we agree that it is necessary for the LSM community to transition to the mixed form of Richards’ equation and to revisit the formulation of the boundary condition.
Line 196: We agree with the reviewer that virtually all LSMs use nonequidistant grids. In Line 196 we state that most LSMs use a staggered grid. In the manuscript we treat the equidistant grid first to further strip the complexity of the solver. We do not stress the implementation of a nonequidistant grid (again, this is a community standard). However, we stress the coordinate transformation as it can be shown that a naïve finite difference discretization for nonequidistant grids leads to the introduction of first order numerical errors (Sundqvist and Veronis, 1970). In practice, we found these errors to be negligible compared to the errors introduced by inadequate grid spacing, which is why we do not discuss this issue any further.
Line 199: Indeed, in order to resolve the gradients in space, thin layers are required. Yet, as the two are dependent on each other (CFL criterion), this enforces the use of a small timestep to ensure numerical stability, or use of an adequate numerical method such as an implicit solver.
Subsection 3.2 “Experiment A”: This whole manuscript is on numerical convergence. Of course, this does not imply that the underlying physical processes are represented correctly. This is not the scope of a convergence study.
Line 273: In the case of a fully saturated soil column, the infiltration is further limited by the outflow at the bottom of the soil column, which is in our case given by the parameterization of Q.
Θ_{f} is relative moisture content defined according to the equation in Table 2.
Subsection 3.3 “Experiment B”: We believe that the description of the forcing is important in this context as it directly relates the numerical errors to different setups that are typical for weather and climate modelling. In our opinion, it is important to show that in a practical context, numerical errors tend to aggravate with precipitation intensity because of the current evolution of weather and climate models towards higher resolutions. The different precipitation time series are chosen, such that they represent a typical example in the transition from convection parameterized to convection resolving models with the associated change in precipitation intensity.
Subsection 4.1 “Analytical Solution” (IA): In the presence of a sink term (Q), it is possible to define a zeroflux bottom boundary condition for K (Schlemmer et al., 2018). In a fully saturated soil, Q is then balanced by infiltration I. This balance controls the magnitude of the flux (i.e. K as D=0) in the soil column. This equilibrium may only be established, if the solver is mass conserving.
References
Ban, N., Schmidli, J., and Schär, C. (2014). Evaluation of the convectionresolving regional climate modeling approach in decadelong simulations. Journal of Geophysical Research: Atmospheres, 119(13):7889–7907.
Balsamo, G., Beljaars, A., Scipal, K., Viterbo, P., van den Hurk, B., Hirschi, M., and Betts, A. K.: A revised hydrology for the ECMWF model: Verification from field site to terrestrial water storage and impact in the Integrated Forecast System, Journal of Hydrometeorology
Decharme, B., Douville, H., Boone, A., Habets, F., and Noilhan, J.: Impact of an exponential profile of saturated hydraulic conductivity within the ISBA LSM: simulations over the Rhône basin, Journal of Hydrometeorology, 7, 61–80, 2006
Ducharne, A.: The hydrol module of ORCHIDEE: Scientific documentation, 2016.
Heise, E., Lange, M., Ritter, B., and Schrodin, R. (2003). Improvement and validation of the multilayer soil model. COSMO newsletter, 3:198–203.
Ogden, F. L., Allen, M. B., Lai, W., Zhu, J., Seo, M., Douglas, C. C., and Talbot, C. A.: The soil moisture velocity equation, Journal of Advances in Modeling Earth Systems, 9, 1473–1487, https://doi.org/https://doi.org/10.1002/2017MS000931, 2017.
Schlemmer, L., Schär, C., Lüthi, D., & Strebel, L. (2018). A groundwater and runoff formulation for weather and climate models. Journal of Advances in Modeling Earth Systems, 10(8), 18091832.
Šimůnek, J., Van Genuchten, M. T., & Šejna, M. (2016). Recent developments and applications of the HYDRUS computer software packages. Vadose Zone Journal, 15(7), vzj201604.
Sundqvist, H. and Veronis, G.: A simple finitedifference grid with nonconstant intervals, Tellus, 22, 26–31, 1970.
Vereecken, H., Weihermüller, L., Assouline, S., Šimůnek, J., Verhoef, A., Herbst, M., Archer, N., Mohanty, B., Montzka, C., Vanderborght, J., Balsamo, G., Bechtold, M., Boone, A., Chadburn, S., Cuntz, M., Decharme, B., Ducharne, A., Ek, M., Garrigues, S., Goergen, K., Ingwersen, J., Kollet, S., Lawrence, D.M., Li, Q., Or, D., Swenson, S., de Vrese, P., Walko, R., Wu, Y. and Xue, Y. (2019), Infiltration from the Pedon to Global Grid Scales: An Overview and Outlook for Land Surface Modeling. Vadose Zone Journal, 18: 153 180191. https://doi.org/10.2136/vzj2018.10.0191

AC2: 'Reply on RC3', Daniel Regenass, 09 Dec 2021

CC1: 'Comment on hess2021426', Ageeth de Haan, 02 Nov 2021

AC4: 'Reply on CC1', Daniel Regenass, 10 Dec 2021
We thank Ms de Haan for her careful and thorough review in the form of a community comment. In the following, we will address some of the main points that are raised in the review. We appreciate that Ms de Haan has found our paper to be overall clearly written and thank her for her positive feedback on that matter.
We are not sure if we understood the first major comment correctly: The interaction of rainfall intensities and numerical errors is one of the main points in this manuscript. However, we fully acknowledge that this should be made clearer in the introduction. Indeed, the main problem with numerical errors of Richards’ equation in the context of weather and climate models is precisely that precipitation intensities increase in highresolution setups. This is a consequence of three main factors: First, increasing horizontal resolution decreases the spatial smoothing of precipitation intensities. Second, shorter time steps are required in order to ensure numerical stability as implied by the CourantFriedrichsLewy (CFL) condition, this in turn decreases the temporal smoothing of precipitation intensities. Third, high resolution setups are often run without a parameterization for convective processes in the atmosphere, which also typically leads to higher precipitation intensities (see e.g. Ban et al., 2014). In that sense the ‘numerical demons’ described in the relatively new paper by la Follette et al. (2021) are also summoned by the recent trend towards kilometerscale modelling in weather and climate models and there are some interesting parallels between the work of la Follette et al. (2021) and our manuscript. As Ms de Haan mentions, the scope of the paper by la Follette et al. (2021) is on conceptual hydrological models while our manuscript is solely on numerical errors of the Richards equation. We choose this relatively narrow scope because the onedimensional Richards’ equation is extensively used in land surface modelling. We acknowledge that the relationship between numerical errors and precipitation intensities could be investigated more systematically.
We feel that the following three paragraphs may be addressed together. There seems to be some misunderstanding concerning the coupling of precipitation on the one hand and Richards’ equation as a part of the land surface model on the other hand. Most modern weather and climate models solve the equations for atmospheric motion on a threedimensional grid (see e.g. Baldauf, 2011), while the land surface model (mainly providing the lower boundary condition for the atmospheric model) is typically onedimensional (see e.g. Balsamo et al., 2009). This means that the rainfall from an atmospheric grid cell with horizontal extent of several kilometers is directly coupled to the onedimensional Richards equation. In that sense, the solution to Richards’ equation may be seen as the grid cell average downward water propagation in the land surface model. Whether or not this is an adequate representation of the underlying physical system is an interesting research question in itself and there are of course good reasons to believe that the Richards equation is not applicable for such large scales (Beven and Cloke, 2012). However, here we solely focus on numerical aspects of Richards’ equation in order to keep the scope of the manuscript as clear as possible. When we refer to numerical issues that arise in kilometerscale modelling, we mainly mean their interaction with high precipitation intensities as such precipitation intensities are simply not present in atmospheric models with a lower resolution. In order for the land surface models to keep up with the trends towards kilometerscale modelling, these issues must be addressed. On a different note, there are approaches to solve the fully threedimensional Richards equation on continental scales (Maxwell et al, 2015). However such approaches require a considerable amount of computing power which is also required for other parts of fully coupled weather and climate models. Mainly due to the computational constraints, the weather and climate community is so far reluctant to move towards a threedimensional representation of water transport in the soil matrix. However, we fully acknowledge that more research is required in order to find suitable approximations for soil water transport on horizontal scales of the order of one kilometer.
References
Baldauf, M., Seifert, A., Förstner, J., Majewski, D., Raschendorfer, M., and Reinhardt, T.: Operational convectivescale numerical weather prediction with the COSMO model: description and sensitivities, Monthly Weather Review, 139, 3887–3905, 2011.
Balsamo, G., Beljaars, A., Scipal, K., Viterbo, P., van den Hurk, B., Hirschi, M., and Betts, A. K.: A revised hydrology for the ECMWF model: Verification from field site to terrestrial water storage and impact in the Integrated Forecast System, Journal of Hydrometeorology
Ban, N., Schmidli, J., and Schär, C. (2014). Evaluation of the convectionresolving regional climate modeling approach in decadelong simulations. Journal of Geophysical Research: Atmospheres, 119(13):7889–7907.
Beven, K. J. and Cloke, H. L.: Comment on “Hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth’s terrestrial water” by Eric F. Wood et al., Water Resources Research, 48, https://doi.org/https://doi.org/10.1029/2011WR010982, 2012.
La Follette, P. T., Teuling, A. J., Addor, N., Clark, M., Jansen, K., & Melsen, L. A. (2021). Numerical daemons of hydrological models are summoned by extreme precipitation. Hydrology and Earth System Sciences, 25(10), 5425–5446. https://doi.org/10.5194/HESS2554252021
Maxwell, R. M., Condon, L. E., and Kollet, S. J.: A highresolution simulation of groundwater and surface water over most of the continental US with the integrated hydrologic model ParFlow v3, Geoscientific Model Development, 8, 923–937, https://doi.org/10.5194/gmd8923 2015, 2015.

AC4: 'Reply on CC1', Daniel Regenass, 10 Dec 2021
Status: closed

RC1: 'Comment on hess2021426', Jan Hopmans, 30 Sep 2021
Initially, when accepting to review this manuscipt, I was very interested to learn about more recent advances in numerical modeling. However, it took me some time to realize that the manuscript assumes knowledge at the landscape scale, and is seeking to describe LS modeling to the km scale, at much higher spatial and temporal resolution. This reviewer's expertise is in soil and vadose zone hydrology which simulates subsurface water flow using the Richards equation at the cm (m) and hour (minute) scales.
I do fully understand the inherent compiexities when taking the landscape scale as the benchmark and using those experiences to address unsaturated water flow at much higher spatial and temporal resolutions. However, over the past few decades, the soil/vadose zone hydrology field has largely advanced, and has developed computer simulation algorithms that can address these high resolutions with satisfactory computer times including multidimensional flow at increasing space (time) scales (more than cm/minute scales).
It was therefore surprising to learn that the authors have largely neglected to reference that work, and ignored those more recent developments of the past 1020 years in the soils and vadose zone hydrology literature. Instead, this manuscript makes simplistic assumptions that are no longer necessary, such as considering multidimensionality through an empirical sink term, and simple boundary conditions such as zero flux lower boundary and a nondynamic infiltration term.
Just for the author’s information, I list hereby just a few references that they could read (from a much larger list of additional references), and these are:
Recent Developments and Applications of the HYDRUS Computer Software Packages
JiÅí ŠimÅ¯nek,Martinus Th. van Genuchten,Miroslav Šejna,at https://acsess.onlinelibrary.wiley.com/doi/full/10.2136/vzj2016.04.0033
Adaptive Grid Refinement in Numerical Models for Water Flow and Chemical Transport in Soil: A Review R. S. Mansell,Liwang Ma,L. R. Ahuja,S. A. Bloom,
https://doi.org/10.2136/vzj2002.2220 at: https://acsess.onlinelibrary.wiley.com/doi/abs/10.2136/vzj2002.2220
Sustainability of irrigated agriculture in the San Joaquin Valley, California
Gerrit Schoups, Jan W. Hopmans, Chuck A. Young, Jasper A. Vrugt, Wesley W. Wallender, Ken K. Tanji, and Sorab Panday at: https://www.pnas.org/content/102/43/15352
And
Modeling Soil Processes: Review, Key Challenges, and New Perspectives
Vereecken,* A. Schnepf, J.W. Hopmans, M. Javaux, D. Or, T. Roose,et al. Modeling Soil Processes: Review, Key Challenges, and New PerspectivesH. Vereecken,* A. Schnepf, J.W. Hopmans, M. Javaux, D. Or, T. Roose et al. 2016. at: https://acsess.onlinelibrary.wiley.com/doi/epdf/10.2136/vzj2015.09.0131

AC1: 'Reply on RC1', Daniel Regenass, 30 Nov 2021
We thank Jan Hopmans for his review and the suggested literature on modelling vadose zone processes. In particular, we read with great interest the article on recent developments in the HYDRUS software package by Šimůnek et al. (2016) as it gives us a good insight on what is considered stateoftheart in the vadose zone community.
We believe that the reviewer (along with reviewer 3) might have misunderstood the aim of this paper and we would greatly appreciate feedback on how to reformulate the manuscript to make its aim more clear. It is not the intent of the authors to suggest a novel numerical method to solve the Richards equation. In particular, the solver presented here has already been introduced by Schlemmer et al. (2018) and it is based on the solver used in the TERRA ML land surface model (LSM) (Heise, 2003) and we hope that this has been made transparent in sections 1 and 2 of our paper. The LSM community and (as part of it) the authors are aware of the progress that has been made in the numerical solution of the Richards equation in the last twenty years. As an example, we explicitly mentioned the recent implementation presented in Lawrence et al. (2019) based on the method introduced by Kavetski et al. (2001), which is comparable to what is used in the HYDRUS model (Šimůnek, 2009). We also referenced publications that address a threedimensional implementation (Maxwell et al. 2015, Mastrotheodoros et al. 2020).
Note that the current version of the CLM model (Lawrence et al., 2019) is a noteworthy exception. In its most recent implementation, the numerical solution of the Richards equation uses an adaptive time step. In contrast, virtually all European centers and modelling consortia use a solver that is similar to the one found in the TERRA ML LSM (Heise, 2003), and use the same time step as in the atmospheric model. This paper is an attempt to raise awareness in the weather and climate modelling community that the numerical methods in the corresponding LSMs are often insufficient to achieve numerical convergence of the Richards equation given high precipitation intensities. This is an urgent matter, as higher precipitation intensities are a natural consequence of the transition in weather and climate modelling to the kilometer scale (see e.g. Ban et al., 2014).
Clearly, LSMs should be adapted to face this challenge and one evident way forward would be to 'catch up' with the progress made in vadose zone physics in the last twenty years. While this is certainly a reasonable and valid approach, we will in the following take the HYDRUS 1D model as an example to highlight potential issues that would arise when using a comparable model as a module inside a weather/ climate model. Eq. 7.1 in Šimůnek et al. (2009) is solved with a Picard iterative solution scheme. Following the notes on discretization (Šimůnek, 2021) this means in practice that up to 10 iterations per time step are required to achieve convergence in the presence of sharp gradients in hydraulic head (which is the case during the infiltration process, in particular during convective events in summer). Moreover, again following the notes on discretization (Šimůnek, 2021), in the presence of sharp gradients, the vertical discretization should be in the order of cm and the 'stretch factor' between two vertical layers should not exceed 1.5. Note that these recommendations are pretty much in line with the findings presented for the Schlemmer et al. (2018) solver investigated in our manuscript in a clear and systematic manner  with the exception that a stretch factor of 1.5 would probably be too high for the scenarios presented, as it would roughly correspond to the grid layout with 11 layers. We are concerned with the application of such a solver in fully coupled weather and climate models. For convection resolving simulations, implementing an iterative solver for the Richards equation would on average probably increase the time to solution of the hydrology module by a factor of three to five at most as the time steps in these models are usually shorter than one minute. Compared with the overall time to solution of the fully coupled model this would probably be acceptable but not optimal as variations in the time to solution are largely unwanted in timecritical applications such as numerical weather prediction, especially in operational contexts. More issues would arise due to the higher vertical resolution as this would not only increase time to solution but also the memory footprint of the model. Given our results and the fact that most soil models have soils deeper than 2 m, it seems unlikely to achieve an acceptable convergence with less than 3040 vertical layers. As highlighted in the manuscript, this is roughly a factor 34 more than what most LSMs use as a standard setup. Given that, it is uncertain, whether the Richards equation captures the relevant physics on the scales of interest (Beven and Cloke, 2012) it is not evident whether the weather/ climate modelling community would be willing to invest that amount of resources on the solution of the Richards equation. In this context, our manuscript could serve the weather/ climate modelling community as a starting point for a discussion on potential ways forward.It is also worth pointing out that the same modeling framework is often used for horizontal resolutions ranging from O(100 km) to O(1 km), increasingly also in global climate models (Stevens et al. 2020). This wide range of resolutions implies significant challenges both on conceptual and computational levels.
References
Ban, N., Schmidli, J., and Schär, C. (2014). Evaluation of the convectionresolving regional climate modeling approach in decadelong simulations. Journal of Geophysical Research: Atmospheres, 119(13):7889–7907.
Beven, K. J. and Cloke, H. L. (2012). Comment on “Hyperresolution global land surfacemodeling: Meeting a grand challenge for monitoring Earth’s terrestrial water” by Eric F. Wood et al. Water Resources Research, 48(1).
Heise, E., Lange, M., Ritter, B., and Schrodin, R. (2003). Improvement and validation of the multilayer soil model. COSMO newsletter, 3:198–203.
Kavetski, D., Binning, P., and Sloan, S. (2001). Adaptive time stepping and error control in a mass conservative numerical solution of the mixed form of Richards equation. Advances in Water Resources, 24(6):595–605.
Lawrence, D. M., Fisher, R. A., Koven, C. D., Oleson, K. W., Swenson, S. C., Bonan, G.,Collier, N., Ghimire, B., van Kampenhout, L., Kennedy, D., et al. (2019). The Community Land Model version 5: Description of new features, benchmarking, and impact of forcing uncertainty. Journal of Advances in Modeling Earth Systems, 11(12):4245–4287.
Schlemmer, L., Schär, C., Lüthi, D., & Strebel, L. (2018). A groundwater and runoff formulation for weather and climate models. Journal of Advances in Modeling Earth Systems, 10(8), 18091832.
Šimůnek, J., Van Genuchten, M. T., & Šejna, M. (2016). Recent developments and applications of the HYDRUS computer software packages. Vadose Zone Journal, 15(7), vzj201604.
Šimůnek, J., Van Genuchten, M. T., & Sejna, M. (2009). The HYDRUS1D software package for simulating the onedimensional movement of water, heat, and multiple solutes in variablysaturated media. University of CaliforniaRiverside Research Reports, Version 4.08.
Šimůnek J. (unknown date). Notes on Spatial and Temporal Discretization (when working with HYDRUS). Retrieved from http://www.pcprogress.com/Documents/Notes_on_Spatial_and_Temporal_Discretization.pdf, on November 14, 2021
Stevens, B. and et al., 2020: The Added Value of Largeeddy and Stormresolving Models for Simulating Clouds and Precipitation. J. Meteorol. Soc. Japan, 98 (2) , 395435 https://doi.org/10.2151/jmsj.2020021

AC1: 'Reply on RC1', Daniel Regenass, 30 Nov 2021

RC2: 'Comment on hess2021426', Anonymous Referee #2, 08 Oct 2021
The comment was uploaded in the form of a supplement: https://hess.copernicus.org/preprints/hess2021426/hess2021426RC2supplement.pdf

AC3: 'Reply on RC2', Daniel Regenass, 09 Dec 2021
We thank the reviewer for their careful review and their comments on our manuscript. In the following, we will discuss their major comments:
 The first major comment addresses the use of the saturationform of the Richards equation (as used in our paper), as opposed to the headbased formulation. We fully agree with the reviewer that the mixedform of Richards’ equation is to be preferred over the saturation form – mainly because of the issues pointed out in the appendix. We are currently working on an implementation of the mixed form and plan to revisit our analysis once the new solver is ready. It is however important to stress that most European weather centers/ weather and climate modelling consortia use the saturation form of Richards’ equation (see manuscript lines 147150). The aim of this paper is not to suggest a novel approach to the numerical solution of the Richards equation, but rather to raise awareness on the (poor) convergence of Richards’ equation as used in weather and climate models and to shed light on the associated errors. Concerning the upper boundary condition (i.e. the formulation of infiltration): We read the corresponding formulation (Eq. 21) in Tubini and Rigon (2021) and think that our numerical implementation is adequate in the saturation form (Eq. 9) as it is a flux boundary condition prescribed by the precipitation flux. Note however, that in the saturation form, an upper limit for this flux must be prescribed as the available pore volume is finite and the flux cannot exceed saturated hydraulic conductivity, because the suction (diffusive) term must vanish at saturation.
 The term ‘successfully’ might indeed be misleading. Here we mean that land surface models using the saturation form of Richards’ equation are used in land surface models and yield physically meaningful results. This does however not imply that the mathematical/ numerical implementation is correct in a rigorous sense. We will try to clarify this issue in the revised version of the manuscript.
References
Tubini, N., & Rigon, R. (2021). Implementing the Water, Heat and Transport model in GEOframe WHETGEO1D v. 1.0: algorithms, informatics, design patterns, open science features, and 1D deployment. Geoscientific Model Development Discussions, 147

AC3: 'Reply on RC2', Daniel Regenass, 09 Dec 2021

RC3: 'Comment on hess2021426', Anonymous Referee #3, 21 Oct 2021
The authors presented in the manuscript a study dealing on ‘improvements’ in temporal and (vertical) spatial resolution in LSM. From the title alone I was attracted as the title suggest novelty and analysis, whereas the manuscript itself does not show any novelty at all. For me it seems y that some basic soil physical principles are not totally clear to the authors. Additionally, recently, Vereecken et al. (2019) published an overview article (Infiltration from the Pedon to Global Grid Scales: An Overview and Outlook for Land Surface Modeling), where from page 8 onwards exactly the problems raised here are shortly discussed. Additionally, there is a link provided (https://www.pcprogress.com/Downloads/Public_Lib_H1D/Using_HYDRUS1D_to_Simulate_Infiltration.pdf) to a full set of simulations using a physical based model showing the effects of different spatial and temporal discretization on infiltration.
Additionally, there are many flaws and sometimes easy problems are made complicated. I worked over the manuscript in detail to about half and made many suggestions and provide comments either in the attached scan or below. After, I decided to have only a quick look at the results as I already made my opinion that the manuscript in its current form is not appropriate for publication.
ABSTRACT
Line 5: If you want to represent runoff and infiltration correctly, even subhourly data are needed, as short but intensive rainfall (e.g., thunderstorms) often cause high runoff as they often reach dry soils.
Last sentence in the abstract: I do not agree, as this is fact irrespectively of the scale models run at (as I also believe that the models commonly used will not change in their algorithms if smaller horizontal resolution will be used). To my understanding the scale (horizontal resolution) only changes the driving forces, especially precipitation.
INTRODUCTION
Line 63: This is exactly the problem in LSM. Here, it was suggested to use a depthdependent profile of saturated hydraulic conductivities. Maybe by this more water can infiltrate but still we will have a wrong model setup generating good results. For sure, hydraulic conductivity will change over the depth profile (which has been shown in many publications and often Ks is high in ploughed layers) but why not first take care that the infiltration process is described correctly?
Line 7580: To disentangle the effect of how Ks will be defined at different scales and what the effects of different forcing are (here precipitation) a simple scaling exercise could have been performed, such as those done for crop models and water balances by Constantin et al. (2019, doi.org/10.1016/j.agrformet.2019.05.013)
Line 82: by some form of the Richards equation. Exactly, this is the problem, as the Richards equation will be solved in the diffusivity form in most LSM and not in the mixed form as most (small) scale soil hydrological models. And this exactly generated the problems as an ‘external’ infiltration concept has to be implemented to solve the infiltration/runoff partitioning (such as GreenAmpt for example). An overview that this holds for most LSM is provided in Tab. 2 of the Vereecken et al. (2019) already mentioned above.
Line 8687. Indeed, different PTF will allow higher or lower infiltration but with respect to the problem analysed here this is not a solution. Would be the same ‘wrong change’ as already done for the vertical dependent Ks.
Line 8892. Indeed, if high resolution maps are used it does not make sense to use class PTFs for e.g. the 12 (or in some LSM only 11) USDA soil classes, as the same output will be generated as we would have for coarser resolution of the soil maps.
Line 9293: Indeed, as from a soil hydrological perspective this has been already analysed over the last 2030 years but seems not to be percolated to the LSM community (or at least part of this community).
Line 100: Adaptive timestepping. At least in the soil hydrological model I use, adaptive time stepping has been implemented since the early 1990s.
Overall, it would be good just to benchmark the LSMs against analytical solutions of infiltration or against soil hydrological models and one would see if they can capture infiltrationrunoff partitioning, infiltration front propagation etc.
THE SYSTEM UNDER CONSIDERATION
This section is somehow textbook knowledge and remind me of what I show in my lectures!
Line 121123: The features you encountered in the experiment are not necessarily connected to soil heterogeneities. It is known that even homogeneous sandy soils will show so called finger flow, which are caused by the hydrophobicity of the sand grains. One a grain will be wetted, the sand turns to be hydrophilic and water can enter the profile, while other regions will stay dry. Additionally, this section is a bit out of scope here and mixes up different processes and problems with those you are working on.
SOIL HYDRAULIC MODEL
Line 159: Here, you state that the two soil hydraulic models need parameterization, which depend on soil texture. This is not fully correct. They depend on the soil pore (distribution), which depend on one hand on soil texture, but also on soil structure. For the conductivity function it is even more complicated as it also depends on pore connectivity.
Line 159. Where did you get the parameters from for the loamy soil? Which PTF or better class PTF was used? Looking at the parameters presented in Tab. 1 I wonder why you showed m for MvG and not n as m is classicaly bound to n by M = 11/n. I also wonder about the alpha value which is 3.6 and therefore too large for classical loamy soils. Finally, why is the saturated hydraulic conductivity Ks different between the two models used as Ks is a material constant of the soil itself?
GROUND RUNOFF
Maybe I got something totally wrong here but isn’t this Dunne overland flow if the soil profile is saturated (often also denoted as saturation excess overland flow)? If yes, this is totally different (in its physical principles) as runoff generated by Hortonion overland flow, which is the precipitation excess with respect to the infiltration capacity?
Finally, why assuming a zeroflux boundary condition at the lower end. To my understanding and neglecting lateral water flow in the saturated zone this will later or earlier always generate Dunne overland flow if the ET is smaller as precipitation.
BOUNDRAY CONDITION
Here we are. In your case, Imax is set to Ks. In reality, this is only the case if infiltration will be only driven by gravity and not by capillary forces. This might happen if the soil is close to saturation but if the soil is much dryer infiltration at its onset can be larger as Ks. This is simple infiltration theory.
NUMERICAL IMPLEMENTATION FOR EQUIDISTANT GRIDS
Line 196: Are you sure that most LSM solve the Richards equation on equidistant discretization. I do not believe so, as to my knowledge most have finer discretization at the top end and increase the spacing with depth. This is also stated in line 216, where it is mentioned that LSM commonly choose grids with finer spacing close to the surface.
Finally, I do not understand why you have to stress on implementing nonequidistant gridding in your model as this is stateoftheart in all other models.
Line 199: There are two points mixed up. The special and temporal resolution.
EXPERIMENT A…
Line 273: …reproduce the analytical solution. What is the analytical solution? Physically, maximum infiltration at saturated conditions without ponding water at the surface equals Ks.
To not understand what you mean by convergence. Numerical convergence does not necessarily mean that the process is represented physically (e.g. shape and wetting front speed)
EXPERIMENT B
Actually, it is not really of interest here how the forcing was calculated by COSMO as any forcing could have been used to show what you want to show. I would either delete or shift to the Annex (if any Annex is allowed in HESS)
FIGURES
Figure 5: Sorry if you define Theta_f as volumetric water content as done in Tab. 1 how can you reach Theta_f of 1 here. This would mean only water and no soil in the systems! I would denote this as lake or ocean!
TABLES
In general, headers for tables should be above the table and not below (see Tab 1 and 2).

AC2: 'Reply on RC3', Daniel Regenass, 09 Dec 2021
We thank the reviewer for taking their time to review our manuscript. First, we would like to stress that the scope of this paper is not on ‘improvements’ but rather to systematically asses numerical issues with the Richards equation as implemented in today’s stateoftheart land surface models (LSMs) coupled to weather and climate prediction models, and to highlight their interaction with increasing precipitation intensities. Since both reviewer 1 and 3 seem to have understood our manuscript as a suggestion for a new solver, we acknowledge that the actual scope has to be made clearer in the introduction. Any feedback on how to make the scope and aim of this paper more understandable would be greatly appreciated. It is clear to the authors that the stateoftheart in many if not most land surface models (e.g. Heise et al., 2003, Balsamo et al., 2008) is vastly different from the stateoftheart in infiltration models used in the vadose zone community (e.g. Šimůnek et al, 2016). Yet, it must be kept in mind that these LSMs are run under operational constraints, with severe limitations on the execution times. This is well documented in the review by Vereecken et al. (2019) and we thank the reviewer for pointing us to this review. While Vereecken et al. (2019) touch on the subject of numerical implementation, they do not discuss numerical convergence analyses, and nor the role of such analyses for establishing the magnitude of respective numerical errors. We thus still believe that systematic convergence analyses of the Richards equation yield added value for the discussion on how to treat the infiltration problem in LSMs. Moreover, the interaction of numerical errors in the solution to Richards equation with precipitation intensities in the context of weather and climate modelling has so far received little attention and are this topic is thus novel for the LSM community. Systematic analyses of numerical errors give us a better baseline for the decision, what to spend the limited execution time on.
In the following, we will address the individual points raised by the reviewer. The reviewer also mentions flaws that they pointed out in handwritten annotations in the manuscript (see supplement to RC3). There is a considerable overlap between the points raised directly in the review and the handwritten annotations in the supplement. We will thus try to address the handwritten comments, where they are not typos/ formal issues and where we feel that they do not correspond to a point raised in the review and mark the corresponding answers with (IA) for ‘inline annotation’.
Line 5: Yes but the feedback to the atmosphere will take place on hourly to seasonal time scales. This is what the weather/ climate modelling community is mostly concerned with.
Line 10 (IA): Here we mean grid spacing.
Line 17 (IA): We think that this will become clear from the experiments, but we could say fast increase in layer thickness in this case means exponential increase in layer thickness as is the case for geometric grids (Table 3).
Line 2123: The solver in the LSM will not change when changing the horizontal resolution of the LSM. But in a coupled setup, the LSM will be subject to a different precipitation forcing. This stems from (a) the effect of smoothing of the simulated precipitation in a coarser atmospheric model and (b) the effect of the parameterization of convection. While parametrized convection in largescale models tends to produce more drizzle and less highintensity precipitation, in kilometerscale models the convective precipitation is simulated explicitly leading to more highintensity precipitation. Both effects are further explained in the introduction and disentangled in chapter 4.3.
Line 4854 (IA): The equations of a dynamical core in the atmosphere are strictly based on conservation laws (mass, momentum, energy) that are described by some simplification of the NavierStokes equations. We do not have an according set of equations to describe conservation laws for the land surface. Therefore, land surface modelling (along with other subgridscale parameterizations in weather and climate models) is subject to approaches that are more heuristic. While this is reasonable given the complexity of the matter, it makes the interactions of the equations at hand (and their numerical representations) harder to understand and quantify.
Line 63: We deliberately choose an implementation of Richards equation that is used in a weather and climate model (Schlemmer et al.,2018). Depth dependent profiles of hydraulic conductivity may be found in similar forms in other LSMs (e.g. Decharme, 2006; Ducharne, 2016). We fully agree that infiltration processes are often poorly represented in weather and climate models. We are however convinced that the solution to the governing equations must be represented accurately by the solver, and to this end convergence studies such as the one presented in our manuscript are important and should be carried out complementary to the implementation of additional processes and a better process representation.
Line 7580: Again, this is mostly a problem of the forcing precipitation data in the coupled model. Here we are specifically concerned with the infiltration process in weather and climate models, where rainfall intensities are not solely subject to spatial aggregation (i.e. the size of the grid cell or LSM tile in case the rainfall is disaggregated), but also the way precipitation is represented. Mesoscale weather models and regional climate models are often operated on the kilometerscale, where convective processes need not to be parameterized in contrast to global models. There are marked differences in precipitation intensities between socalled convection permitting (or convection resolving) models and models with parameterized convection processes (Ban et al., 2014). The scaling of hydraulic properties such as hydraulic conductivity is beyond the scope of this paper, but we agree that the subgridscale heterogeneity should urgently be addressed in future research. Thank you for pointing us to the reference.
Line 82: We fully agree with the reviewer that the mixed form of Richards equation should be used and we are currently revisiting our analysis for the mixed form. However, as can be seen in Vereecken et al. (2019), Table 2, four out of seven of the listed LSMs use the diffusivity (saturation) form of Richards’ equation. Moreover, five out of seven of the listed LSMs generate surface runoff when the precipitation rate exceeds saturated hydraulic conductivity. Therefore, a convergence analysis on such an implementation is of practical relevance to the LSM community.
Lines 8687 and 8892: There are ongoing attempts to improve the representation of hydrological processes along several different pathways. None of those will be efficient on its own. We agree with the reviewer that the infiltration process is a particularly critical one, and it must first be represented accurately by the corresponding governing equations and their numerical representation. Otherwise errors in the dataset will interact with physical and numerical errors, which hampers LSM development.
Lines 9293: The lack of awareness in the LSM community – or, more broadly speaking, the weather and climate community – is precisely the issue we want to address with our manuscript. This would allow for a discussion on ways forward. We agree with the reviewer that this ‘knowledge gap’ should be closed in order to foster the collaboration amongst hydrologists, soil physicists and LSM modelers.
Line 100: Variations in timetosolution is problematic in the operational context of weather and climate modelling. This might at least partially be the reason, why parts of the LSM community are hesitant to adopt such approaches.
Line 133 and line 150 (IA): In TERRA ML (Heise et al., 2004, Schlemmer et al., 2018), E is a sink term in the root zone. We agree with the reviewer that it would be more elegant to include transpiration in the suction head (if this is what they mean). However most LSMs are simply not that far in their evolution yet.
Line 146 (IA): Q is a sink term that should emulate horizontal transport in complex terrain to some extent. We follow the implementation of Schlemmer et al. (2018). While this approach is of course very heuristic, it keeps the LSM onedimensional and it is computationally efficient.
Benchmarking LSMs against infiltration models is an excellent idea and we are very open to this suggestion.
Subsection 2.2.1 “The system under consideration”: We are open to omit this section. We believe however that it nicely illustrates the problem with vertical resolution. If one compares the sharp moisture gradient across the wetting front with the vertical discretization in stateoftheart LSMs, the issue as well as the effect on infiltration should be immediately clear.
Subsection 2.2.2 “Soil hydraulic model” (line 159): The values are adopted from the ORCHIDEE LSM (Ducharne, 2016). We are open to revisting the analysis with different values, but it will likely not change the results of the convergence analysis tremendously.
Subsection 2.2.3 “Ground Runoff”: The parameterization of ground runoff in our version of TERRA ML (Heise et al., 2003) is described in Schlemmer et al. (2018). In contrast, both saturation excess and infiltration excess are included in the surface runoff. In our case it would correspond to a saturation of the upper soil and an exceedance of saturated hydraulic conductivity respectively.
Subsection 2.2.4 “Boundary condition”: The answer to this point is partially included in the answer to the point raised concerning using the mixed form of Richards’ equation (line 82). As pointed out above, five out of seven of the listed LSMs in Table 2 of the review by Vereecken et al. (2019) generate surface runoff when the precipitation rate exceeds saturated hydraulic conductivity. While this approach is likely not valid for very dry soils, it is still of practical relevance as this assumption is included in many if not most LSMs and it is probably also a fairly good approximation for wet and moderately dry soils (see also Ogden et al., 2017). However, we agree that it is necessary for the LSM community to transition to the mixed form of Richards’ equation and to revisit the formulation of the boundary condition.
Line 196: We agree with the reviewer that virtually all LSMs use nonequidistant grids. In Line 196 we state that most LSMs use a staggered grid. In the manuscript we treat the equidistant grid first to further strip the complexity of the solver. We do not stress the implementation of a nonequidistant grid (again, this is a community standard). However, we stress the coordinate transformation as it can be shown that a naïve finite difference discretization for nonequidistant grids leads to the introduction of first order numerical errors (Sundqvist and Veronis, 1970). In practice, we found these errors to be negligible compared to the errors introduced by inadequate grid spacing, which is why we do not discuss this issue any further.
Line 199: Indeed, in order to resolve the gradients in space, thin layers are required. Yet, as the two are dependent on each other (CFL criterion), this enforces the use of a small timestep to ensure numerical stability, or use of an adequate numerical method such as an implicit solver.
Subsection 3.2 “Experiment A”: This whole manuscript is on numerical convergence. Of course, this does not imply that the underlying physical processes are represented correctly. This is not the scope of a convergence study.
Line 273: In the case of a fully saturated soil column, the infiltration is further limited by the outflow at the bottom of the soil column, which is in our case given by the parameterization of Q.
Θ_{f} is relative moisture content defined according to the equation in Table 2.
Subsection 3.3 “Experiment B”: We believe that the description of the forcing is important in this context as it directly relates the numerical errors to different setups that are typical for weather and climate modelling. In our opinion, it is important to show that in a practical context, numerical errors tend to aggravate with precipitation intensity because of the current evolution of weather and climate models towards higher resolutions. The different precipitation time series are chosen, such that they represent a typical example in the transition from convection parameterized to convection resolving models with the associated change in precipitation intensity.
Subsection 4.1 “Analytical Solution” (IA): In the presence of a sink term (Q), it is possible to define a zeroflux bottom boundary condition for K (Schlemmer et al., 2018). In a fully saturated soil, Q is then balanced by infiltration I. This balance controls the magnitude of the flux (i.e. K as D=0) in the soil column. This equilibrium may only be established, if the solver is mass conserving.
References
Ban, N., Schmidli, J., and Schär, C. (2014). Evaluation of the convectionresolving regional climate modeling approach in decadelong simulations. Journal of Geophysical Research: Atmospheres, 119(13):7889–7907.
Balsamo, G., Beljaars, A., Scipal, K., Viterbo, P., van den Hurk, B., Hirschi, M., and Betts, A. K.: A revised hydrology for the ECMWF model: Verification from field site to terrestrial water storage and impact in the Integrated Forecast System, Journal of Hydrometeorology
Decharme, B., Douville, H., Boone, A., Habets, F., and Noilhan, J.: Impact of an exponential profile of saturated hydraulic conductivity within the ISBA LSM: simulations over the Rhône basin, Journal of Hydrometeorology, 7, 61–80, 2006
Ducharne, A.: The hydrol module of ORCHIDEE: Scientific documentation, 2016.
Heise, E., Lange, M., Ritter, B., and Schrodin, R. (2003). Improvement and validation of the multilayer soil model. COSMO newsletter, 3:198–203.
Ogden, F. L., Allen, M. B., Lai, W., Zhu, J., Seo, M., Douglas, C. C., and Talbot, C. A.: The soil moisture velocity equation, Journal of Advances in Modeling Earth Systems, 9, 1473–1487, https://doi.org/https://doi.org/10.1002/2017MS000931, 2017.
Schlemmer, L., Schär, C., Lüthi, D., & Strebel, L. (2018). A groundwater and runoff formulation for weather and climate models. Journal of Advances in Modeling Earth Systems, 10(8), 18091832.
Šimůnek, J., Van Genuchten, M. T., & Šejna, M. (2016). Recent developments and applications of the HYDRUS computer software packages. Vadose Zone Journal, 15(7), vzj201604.
Sundqvist, H. and Veronis, G.: A simple finitedifference grid with nonconstant intervals, Tellus, 22, 26–31, 1970.
Vereecken, H., Weihermüller, L., Assouline, S., Šimůnek, J., Verhoef, A., Herbst, M., Archer, N., Mohanty, B., Montzka, C., Vanderborght, J., Balsamo, G., Bechtold, M., Boone, A., Chadburn, S., Cuntz, M., Decharme, B., Ducharne, A., Ek, M., Garrigues, S., Goergen, K., Ingwersen, J., Kollet, S., Lawrence, D.M., Li, Q., Or, D., Swenson, S., de Vrese, P., Walko, R., Wu, Y. and Xue, Y. (2019), Infiltration from the Pedon to Global Grid Scales: An Overview and Outlook for Land Surface Modeling. Vadose Zone Journal, 18: 153 180191. https://doi.org/10.2136/vzj2018.10.0191

AC2: 'Reply on RC3', Daniel Regenass, 09 Dec 2021

CC1: 'Comment on hess2021426', Ageeth de Haan, 02 Nov 2021

AC4: 'Reply on CC1', Daniel Regenass, 10 Dec 2021
We thank Ms de Haan for her careful and thorough review in the form of a community comment. In the following, we will address some of the main points that are raised in the review. We appreciate that Ms de Haan has found our paper to be overall clearly written and thank her for her positive feedback on that matter.
We are not sure if we understood the first major comment correctly: The interaction of rainfall intensities and numerical errors is one of the main points in this manuscript. However, we fully acknowledge that this should be made clearer in the introduction. Indeed, the main problem with numerical errors of Richards’ equation in the context of weather and climate models is precisely that precipitation intensities increase in highresolution setups. This is a consequence of three main factors: First, increasing horizontal resolution decreases the spatial smoothing of precipitation intensities. Second, shorter time steps are required in order to ensure numerical stability as implied by the CourantFriedrichsLewy (CFL) condition, this in turn decreases the temporal smoothing of precipitation intensities. Third, high resolution setups are often run without a parameterization for convective processes in the atmosphere, which also typically leads to higher precipitation intensities (see e.g. Ban et al., 2014). In that sense the ‘numerical demons’ described in the relatively new paper by la Follette et al. (2021) are also summoned by the recent trend towards kilometerscale modelling in weather and climate models and there are some interesting parallels between the work of la Follette et al. (2021) and our manuscript. As Ms de Haan mentions, the scope of the paper by la Follette et al. (2021) is on conceptual hydrological models while our manuscript is solely on numerical errors of the Richards equation. We choose this relatively narrow scope because the onedimensional Richards’ equation is extensively used in land surface modelling. We acknowledge that the relationship between numerical errors and precipitation intensities could be investigated more systematically.
We feel that the following three paragraphs may be addressed together. There seems to be some misunderstanding concerning the coupling of precipitation on the one hand and Richards’ equation as a part of the land surface model on the other hand. Most modern weather and climate models solve the equations for atmospheric motion on a threedimensional grid (see e.g. Baldauf, 2011), while the land surface model (mainly providing the lower boundary condition for the atmospheric model) is typically onedimensional (see e.g. Balsamo et al., 2009). This means that the rainfall from an atmospheric grid cell with horizontal extent of several kilometers is directly coupled to the onedimensional Richards equation. In that sense, the solution to Richards’ equation may be seen as the grid cell average downward water propagation in the land surface model. Whether or not this is an adequate representation of the underlying physical system is an interesting research question in itself and there are of course good reasons to believe that the Richards equation is not applicable for such large scales (Beven and Cloke, 2012). However, here we solely focus on numerical aspects of Richards’ equation in order to keep the scope of the manuscript as clear as possible. When we refer to numerical issues that arise in kilometerscale modelling, we mainly mean their interaction with high precipitation intensities as such precipitation intensities are simply not present in atmospheric models with a lower resolution. In order for the land surface models to keep up with the trends towards kilometerscale modelling, these issues must be addressed. On a different note, there are approaches to solve the fully threedimensional Richards equation on continental scales (Maxwell et al, 2015). However such approaches require a considerable amount of computing power which is also required for other parts of fully coupled weather and climate models. Mainly due to the computational constraints, the weather and climate community is so far reluctant to move towards a threedimensional representation of water transport in the soil matrix. However, we fully acknowledge that more research is required in order to find suitable approximations for soil water transport on horizontal scales of the order of one kilometer.
References
Baldauf, M., Seifert, A., Förstner, J., Majewski, D., Raschendorfer, M., and Reinhardt, T.: Operational convectivescale numerical weather prediction with the COSMO model: description and sensitivities, Monthly Weather Review, 139, 3887–3905, 2011.
Balsamo, G., Beljaars, A., Scipal, K., Viterbo, P., van den Hurk, B., Hirschi, M., and Betts, A. K.: A revised hydrology for the ECMWF model: Verification from field site to terrestrial water storage and impact in the Integrated Forecast System, Journal of Hydrometeorology
Ban, N., Schmidli, J., and Schär, C. (2014). Evaluation of the convectionresolving regional climate modeling approach in decadelong simulations. Journal of Geophysical Research: Atmospheres, 119(13):7889–7907.
Beven, K. J. and Cloke, H. L.: Comment on “Hyperresolution global land surface modeling: Meeting a grand challenge for monitoring Earth’s terrestrial water” by Eric F. Wood et al., Water Resources Research, 48, https://doi.org/https://doi.org/10.1029/2011WR010982, 2012.
La Follette, P. T., Teuling, A. J., Addor, N., Clark, M., Jansen, K., & Melsen, L. A. (2021). Numerical daemons of hydrological models are summoned by extreme precipitation. Hydrology and Earth System Sciences, 25(10), 5425–5446. https://doi.org/10.5194/HESS2554252021
Maxwell, R. M., Condon, L. E., and Kollet, S. J.: A highresolution simulation of groundwater and surface water over most of the continental US with the integrated hydrologic model ParFlow v3, Geoscientific Model Development, 8, 923–937, https://doi.org/10.5194/gmd8923 2015, 2015.

AC4: 'Reply on CC1', Daniel Regenass, 10 Dec 2021
Daniel Regenass et al.
Daniel Regenass et al.
Viewed
HTML  XML  Total  BibTeX  EndNote  

877  323  27  1,227  18  17 
 HTML: 877
 PDF: 323
 XML: 27
 Total: 1,227
 BibTeX: 18
 EndNote: 17
Viewed (geographical distribution)
Country  #  Views  % 

Total:  0 
HTML:  0 
PDF:  0 
XML:  0 
 1