It rains and then? Numerical challenges with the 1D Richards equation in kilometer-resolution land surface modelling

Over the last decade kilometer-scale weather predictions and climate projections have become established. Thereby both the representation of atmospheric processes, as well as land-surface processes need adaptions to the higher-resolution. 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 qual5 ity 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 kilometer-scale, convection-permitting models. First of all, the precipitation-intensity distribution changes to more intense events. Second, the time-step of the numerical integration becomes smaller. The aim of this study is to investigate the numerical convergence of the one-dimensional Richards Equation with respect to the soil hydraulic model, vertical layer thickness and time-step during the infiltration process. Both regular 10 and non-regular (unequally spaced) grids typical in land surface modelling are considered, using a conventional semi-implicit 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 5cm and for time steps greater than 40s, irrespective of the soil hydraulic model. The velocity of the wetting front decreases systematically with increasing time step and decreasing vertical resolution. For non-regular grids, a new discretization based on a coordinate transform is introduced. In contrast to 15 simpler vertical discretizations, it is able to represent the solution second-order accurate. The results for non-regular 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 infiltration-runoff partitioning for precipitation time series charac20 teristic 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 high-resolution counterparts. Our analysis indicates that the problem is particularly pronounced for kilometer-resolution models. 1 https://doi.org/10.5194/hess-2021-426 Preprint. Discussion started: 8 September 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
Kilometer-scale numerical weather prediction (NWP) and climate models have approached operational use on continental to 25 global scales (Prein et al., 2015;Fuhrer et al., 2018;Leutwyler et al., 2017;Stevens et al., 2019;Schär et al., 2020;Wedi et al., 2020). This marks a major breakthrough in atmospheric modelling, as a large part of vertical energy redistribution can be resolved explicitly and it is reasonable to assume that global kilometer-scale modelling could fundamentally enhance our understanding of the climate system (Palmer, 2014). Numerous publications point out the improved representation of precipitation over land in convection resolving models (CRM hereafter) (e.g Ban et al., 2014;Leutwyler et al., 2017;Pichelli et al., 2021;30 Vergara- Temprado et al., 2021), in particular the improved representation of the diurnal cycle and rainfall intensities -a longstanding issue in global convection-parameterized NWP and climate models (CPM hereafter) (Stephens et al., 2010). On the one hand, the extensive computational challenges of continental and global convection resolving modelling are evident (Fuhrer et al., 2014;Schär et al., 2020). On the other hand, it is to be expected that higher resolution setups exhibit higher velocities, mass and energy fluxes throughout the whole modelling system, simply as a consequence of reduced spatio-temporal averag-35 ing. While this yields exciting possibilities for the investigation of extreme events in a changing climate (Vergara-Temprado et al., 2021, e.g.), it challenges established approaches in numerics and subgridscale parameterizations inherited from CPMs (Yano et al., 2018) -such as the land surface model.
The recent advances in atmospheric modelling are complemented by similar tendencies in global land surface modelling. 40 Wood et al. (2011) view the development of global, comprehensive land surface simulations on the order of 1 km as a "grand challenge" for the hydrologic community. Indeed, recent sub-continental to continental scale three dimensional land surface simulations yield realistic features of the terrestrial water cycle (Maxwell et al., 2015;Mastrotheodoros et al., 2020). However, as is the case for atmospheric science (Yano et al., 2018), there is a fundamental lack of scientific understanding of subgrid-scale processes that ought to be parameterized even on the spatial scales envisaged by Wood et al. (2011) in their grand challenge 45 paper, as has also been pointed out by Beven and Cloke (2012).
In contrast to models of the atmospheric dynamics, land-surface models (LSMs) often do not have a closed system of governing equations, and advances in process understanding are often translated into heuristic implementations. This may deem it impossible for the individual modeler to have an overview on the full LSM and it leads to a large amount of free parameter 50 which are not well constrained by theory (Fisher and Koven, 2020). Moreover, we argue here that it may also lead to a not fully transparent coupling of individual subprocesses within the LSM and with the atmosphere. In addition, numerical aspects of LSMs are receiving much less attention and scrutiny than in atmospheric models, presumably because of the important role of semi-empirical concepts. 55 Indeed, this study is partially motivated by challenges the authors faced while trying to test and calibrate a new groundwater model for the TERRA ML LSM (Heise et al., 2003) presented by Schlemmer et al. (2018). When comparing a CPM to a CRM and surface runoff sums as obtained from simulations carried out with the COSMO model (Baldauf et al., 2011) following the same setup as used by Zeman et al. (2021) for May 29 2018, a day characterized by convective activity across the entire European continent. Overall, the CRM simulation exhibits higher precipitation sums than the CPM simulation and, more intricate, overall higher amounts of surface runoff. Just by comparing differences in precipitation to differences in surface runoff (Fig. 1, bottom row), it becomes clear that the relation between the two is not trivial. Differences in surface runoff are more patchy, 70 likely an indication for the important role high-intensity precipitation events and surface heterogeneity play for the observed differences in partitioning of water at the surface. The difference between the two quantities in turn represent the water that infiltrates into the soil, which may influence the storage and thus long-term memory of the soil.
The understanding of relevant physical processes and relevant model parameters is further obscured, as in some LSMs, nu-75 merical and physical approximations of the real-world problem are not fully separated. On the physical side, the challenge lies in the accurate estimation of saturated hydraulic conductivity -or even broader, finding a suitable approach to the infiltration problem on the kilometer-scale. This is then followed by the numerical representation of the governing equation for water transport in the soil matrix (in this case the Richards (1931) equation). This has become evident as the authors observed time step and resolution dependencies in the generation of surface runoff in standalone applications of TERRA ML.

80
As is the case for TERRA ML, most 2nd and 3rd generation LSMs describe water transport in the soil matrix by some form of the one-dimensional Richards (1931) equation (Pitman, 2003). This 'backbone' of the hydrology section has not changed dramatically in current state-of-the art LSMs (Balsamo et al., 2009;Niu et al., 2011;Ducharne, 2016;Lawrence et al., 2019;Reick et al., 2021). Extensions to LSM hydrology have primarily focused on groundwater parameterizations (Niu et al., 2007;85 Schlemmer et al., 2018) and infiltration-runoff partitioning (Liang et al., 1994;Niu et al., 2005). Recently,some effort is invested into incorporating new pedotransfer functions into LSMs and it is clear that the choice of pedotransfer functions affects the model properties substantially (Weihermüller et al., 2021). Modern methods to estimate pedotransfer functions (Schaap et al., 2001;Saxton and Rawls, 2006;Van Looy et al., 2017) have a continuous dependency on soil texture instead of being fixed for discrete soil classes and account for organic matter in the soil. Such pedotransfer functions facilitate the use of high 90 resolution maps of soil properties (most notably SoilGrids250m, Hengl et al., 2017).
In contrast to the efforts in adding more realism to soil physical aspects, numerical aspects of the Richards equation have  received much less attention. Nevertheless, there seems to be some awareness of numerical issues in the LSM community. Already quite early after the implementation of Richards equation into LSMs, Lee and Abriola (1999) pointed out that vertically integrated solutions of the Richards equation -as used by many LSMs at that time and still today -tend to overpredict soil 95 moisture in the uppermost layers and underpredict discharge to the aquifer. Campoy et al. (2013) drastically increased the number of vertical layers to conduct experiments concerning the choice of the bottom boundary condition in ORCHIDEE. Zeng and Decker (2009) tried to reduce truncation errors of the spatial discretization by subtracting the equilibrium soil moisture from the equation. An approach that was later discarded by Lawrence et al. (2019) due to problems pointed out by De Rooij (2010). Lawrence et al. (2019) recently introduced an adaptive time-stepping procedure developed by Kavetski et al. (2001) in order 100 to improve numerical stability of solutions to the Richards equation. Mueller-Quintino et al. (2016) compared three LSMs and found spurious oscillations at low vertical resolutions. However, there seems to be a lack of systematic assessment concerning the magnitude of numerical errors introduced by solving the Richards equation with standard approaches and setups.
In order to elucidate the problem of numerical errors, this paper addresses the following research questions: Do numerical representations of the solution to the 1D Richards Equation as used in land surface models converge? How good is the con-105 vergence? What are the implications for the infiltration process? How do numerical errors interact with rainfall intensities as obtained from CPM and CRM simulations, respectively? What are the implications for kilometer-scale weather and climate modeling?
The remainder of this paper is organized as follows: The model and its numerical representation are introduced in section 2.
The convergence experiments are described in section 3 and the corresponding results are presented in section 4. Finally, some 110 final thoughts and our conclusions are formulated in section 5.

Model
In the first part of this section, the derivation of the one-dimensional Richards Equation is briefly revisited and two different soil hydraulic models are discussed. In the second part, the numerical implementation is presented and a novel coordinate transform for the one-dimensional Richards Equation on non-equidistant grids is introduced. The soil model described in this 115 section is implemented in Python.

The System Under Consideration
In the following sections, a model of one-dimensional water transport in a homogeneous, porous medium will be introduced.
Let us first take a step back in order to get a feeling for the system described by the one-dimensional Richards equation. A related physical experiment is shown in Fig. 2 (although the transport here is not strictly one dimensional, but sufficiently close).

120
There are three main observations to be highlighted: First, even in a setup as homogeneous as depicted in Fig. 2 preferential flow paths seem to form, where water transport happens faster. On the order of LSM grid cells or tiles, the inhomogeneities are substantially larger and should be accounted for, but this is beyond the scope of this paper. Second -and this is important for the numerical solution of Richards equation -there is a sharp gradient in soil moisture from the wet to the dry zone, the so-called Figure 2. Do-it-yourself infiltration experiment with a sand column. Water is poured at the top of the column, and the panels show the infiltration after five minutes (left) and one hour (right). Note the sharp gradient in the transition from wet to dry sand -the so-called wetting front. The propagation velocity of the wetting front determines the infiltration rate and hence the amount of ponding water on the surface. wetting front, which is a characteristic feature of water propagation in porous media. This sharp gradient must be mirrored in 125 some way in the soil hydraulic model and ideally be captured by the spatial and temporal discretization. Finally, notice the water ponding on top of the soil column which almost completely vanishes after one hour of infiltration. The propagation of the wetting front determines the amount of water that may infiltrate the soil.

130
Let us consider an arbitrary element of soil with volumetric soil moisture θ in a Cartesian coordinate system. Mass conservation where F is the volumetric soil water flux in m/s and E is a sink term due to local evapotranspiration. Infiltration is given as a boundary condition for the vertical flux at the surface, i.e. F z (z = 0) = I. Let us now assume that the flux F follows Darcy's with hydraulic conductivity K and hydraulic head h = ψ + z t − z. The hydraulic head in unsaturated soils is the sum of suction head ψ and gravitational head z t − z, where z t is the topographic height and z is the depth beneath the local surface. We introduce the horizontal Nabla ∇ h = (∂/∂x, ∂/∂y) and separate lateral and vertical flow (thereby moving to a downward pointing vertical coordinate z) to obtain Due to computational constraints, lateral flow is often neglected or parameterized. Here, lateral flow is treated as a source/sink term Q = −∇ h · K∇ h (ψ + z t − z) and we will discuss its parameterization in section 2.2.3. Using this simplification, applying the chain rule and introducing hydraulic diffusivity D(θ) = K(θ)∂h/∂θ, we can rewrite (3) as (Hillel, 2012) 145 This is the one-dimensional saturation-based form of Richards Equation which shall be used subsequently. As a number of LSMs in operational NWP and climate models use the saturation form of the Richards equation (Heise et al., 2003;Balsamo et al., 2009;Ducharne, 2016;Reick et al., 2021, mainly at European weather services and research institutions, see), we will base our analysis on it. In land surface modelling, evapotranspiration is usually estimated in a different part of the model which 150 is beyond the scope of this paper. Therefore, we set E = 0 for all further derivations and experiments in this study. Still missing are expressions for K, D and Q, which we will introduce in the following subsections.

Soil Hydraulic Model K, D
Following Shao and Irannejad (1999), the functions K and D will be referred to as the soil hydraulic model. Here, two different soil hydraulic models are tested. On the one hand, the Rijtema (1969) model (RT model hereafter) with exponential 155 dependencies on saturation which is used in TERRA ML (Heise et al., 2003), and on the other hand the Mualem (1976) -Van Genuchten (1980) model (MVG model hereafter) with stronger dependencies on soil moisture close to saturation and air dry point, which is preferred by many soil physicists (Shao and Irannejad, 1999). Both models require a number of parameters which depend on soil texture. In this study, parameters are chosen for a loam soil. All parameters and the hydraulic functions are summarized in Table 1. There is one important detail in the MVG model which requires some consideration: The hydraulic 160 diffusivity D diverges as the relative moisture θ f approaches saturation (see A and Van Genuchten (1980)). In practice, D has to be limited to finite values. In ORCHIDEE, this is achieved by piecewise discretization of the function D(θ) with respect to θ. The value of the last interval (with the relative moisture content θ f = 1 as upper bound) is then set to a fraction of the second but last interval (Ducharne, 2016). In HTESSEL (Balsamo et al., 2009) and JSBACH (Reick et al., 2021), D is modeled after the equation given by Clapp and Hornberger (1978), while K is still given by the MVG model. This is somewhat contradictory

Ground Runoff Q
In order to parameterize the effect of topography on ground runoff, Schlemmer et al. (2018) implemented a formulation for ground runoff from the saturated part of the soil which should in principle be independent from the horizontal resolution of the land surface model. They assume runoff divergence (the horizontal contribution to (3)) from the saturated part of the soil 175 column (from bedrock z B to the upper limit of the water table z wt ) to be given by where S ORO is a slope parameter representing the sub-gridscale gradient of orography and L −1 g is the inverse of the dominant length scale of sub-gridscale heterogeneity. Note that the part K 0 · S ORO is equivalent to Darcy flow from the saturated zone.
In this formulation Q is strictly positive due to the specific formulation of S ORO (see Schlemmer et al. (2018) for details).

180
The rationale of this approach is to include the water table in the soil column and to diagnose the position of the water table and subsequently ground runoff. The water table has to be diagnosed each time step, the procedure takes water from the saturated zone as well as from the layer above into consideration, the details are given in Schlemmer et al. (2018). Note that this parameterization requires a zero-flux boundary condition in order for the water table to build up.

Boundary Conditions
The upper boundary condition is given by infiltration I, which is on the one hand given by downward water flux at the surface (precipitation, snow melt, etc.) denoted with P hereafter. On the other hand I is controlled by maximum infiltration I M AX 190 which cannot exceed saturated hydraulic conductivity at the surface I M AX ≤ K 0 (z = 0) and furthermore it is limited by available pore space at the surface. The former limitation is straightforward to implement, the latter is subject to the numerical implementation. Surface runoff is then given by the difference of downward water flux at the surface and infiltration,

195
As is the case for many land surface models (see e.g. Ducharne, 2016;Lawrence et al., 2019;Heise et al., 2003), Richards Equation is solved on a staggered grid with the fluxes on the layer interfaces (Fig. 3). Here, the layers are chosen to be equidistant, a method for non-equidistant grids will be introduced in section 2.4. In order to capture the sharp gradients in soil moisture around the wetting front, relatively thin layers are required. Therfore, to ensure numerical stability, a semi-implicit time stepping procedure is chosen. The diffusive part of (4) is solved with a simplified backward-implicit time stepping. The simplification hereby is that the functions K and D use the volumetric soil moisture of time level m (instead of m + 1). Using the presented spatial discretization and the semi-implicit time-stepping, the discretized form of (4) reads Solving for θ m+1 k on all levels k leads to a tridiagonal linear system of equations which can be solved using standard methods.
Note that the explicit part of (8) uses flux corrected quantities for K and Q as introduced by Schlemmer et al. (2018) following the flux-corrected transport (FCT) method after Boris and Book (1973) and Zalesak (1979). This ensures mass conservation even in the fully saturated part of the soil and close to the air dry point. A special case of the flux correction for overfill (see whereK n 1/2 andQ n 1 are (flux corrected) gravitational flux across the bottom of the first layer and ground runoff from the first layer respectively and I m K = min (P, K 0 (z = 0)) is conductivity limited infiltration (see section 2.2.4).

215
In land surface modelling, it is common to choose a grid with thinner layers close to the surface to capture the infiltration more accurately. At the same time the coarser grid at larger depths allows deeper soils without becoming computationally too expensive. While this seems reasonable at first glance, it should be stressed that (8) (Sundqvist and Veronis, 1970). In this section, (8) is transformed to allow for geometric grids while retaining second-order accuracy 220 ('geometric' as the layer thicknesses can be expressed as a geometric series).

Transformed Equation
Let us introduce ζ such that z = f (ζ), where f is continuously differentiable on the region of interest and the inverse function is continuously differentiable. Using the chain rule, the derivatives in the governing equation (5) can be expanded to yield

225
The transformed equation (10) can not be rearranged into flux form and it would thus be incompatible with the FCT method.
To overcome this issue, (10) is multiplied with the inverse of the metric term J −1 = ∂z/∂ζ to obtain, where J = ∂ζ/∂z is the metric term and the 'ˆ' indicates transformed quantities, e.g.θ = J −1 θ. For the implementation, it is important to note that θ appears in its transformed and non-transformed state and some care has to be taken to avoid 230 inconsistencies when setting up the linear system of equations.

Suggestion for a Geometric Grid
The aim of the transformation is to put more layers close to the surface and fewer layers in the deep soil, where the time scales are often assumed to be larger and the spatial gradients of soil moisture less sharp. An exponential increase of layer thickness seems to be a reasonable approach and we use This transformation has two parameters, the dimensionless constant b that controls the increase in thickness from one layer to the next and ∆z 1 , the thickness of the uppermost layer. Upon requesting ζ (∆z 1 ) = ∆ζ using 12, the equation connecting b and ∆ζ can be derived as It is then straightforward to show that the layer thickness of layer k is given by and using the summation rule for geometric series, the depth of nth layer interface is located at For the limit b → 0, the regular grid is recovered if ∆ζ is chosen according to (13). To complete the transformation rules for 245 this specific transformation, we calculate the metric term at layer k to be and its inverse The metric terms for the layer interfaces may be calculated analogously.

Description of the Numerical Experiments
All experiments are carried out with a soil depth of 2 m meaning that the first layer interface (i.e. the surface) is at z = 0 m and the last interface is at z = 2 m. An initially moist environment is chosen with θ init = 0.8η. Parameter choices related to the soil hydraulic model are given in Table 1. Evapotranspiration is switched off for all experiments in order to focus solely on infiltration and drainage processes.  Table 3.  Table 3. Grid layouts for the experiments with geometric grids. nz is the number of vertical layers, b is the stretch factor and ∆z indicates the range of layer thickness.

Experiment A: Convergence Analysis in an Idealized Infiltration Setup
In the first experiment, the soil column is forced with a constant precipitation rate of P = 20 mm h −1 . This rate exceeds the 265 saturated hydraulic conductivity for both the MVG and the RT soil hydraulic model. Thus, the upper soil saturates quickly, surface runoff sets in, and a sharp wetting front forms subsequently in the soil. Table 2 gives an overview on the model setups tested. The precipitation forcing is kept up until the soil column is fully saturated for all tested setups and then run with the same forcing up to a total of 33.6 hours. The forcing is then stopped and the simulations are subject to a drainage processes steered entirely by groundwater runoff (see section 2.2.3) for the same amount of time, leading to a total model runtime of 67.2 270 hours. Experiment A run on a regular grid is denoted AR, while on a geometric grid it is marked AG.
This experiment has two different purposes: The first one is to test whether the model can accurately reproduce the analytical solution for a fully saturated soil column. This allows for checking whether the model is running correctly. The second purpose is to check, whether the numerical approximation of the infiltration process converges, that is to test whether the propagation 275 of the wetting front depends on the choice of the numerical grid and the time step. If the solution converges, it is possible to define a range of acceptable choices for grid setup and time step respectively.
The simulation with 200 equidistant soil layers of ∆z = 1 cm thickness and with a temporal resolution of ∆t = 1 s is defined as the reference simulation and numerical convergence is always investigated with respect to this reference. As a 280 metric of convergence, normalized root mean squared differences of surface runoff Q S with respect to the reference simulation are calculated by

Experiment B: Convergence Under Realistic Precipitation Forcing
As experiment A is subject to a rather unrealistic forcing, a second experiment with a more realistic precipitation forcing is km simulation and explicitly considered in the 2.2 km simulation. The shallow convection parameterization scheme (Tiedtke, 1989) is switched on for both simulations. Fig. 1 shows the 24 h sum of precipitation as obtained from the two simulations. One local maximum which is present in both simulations is located at the Southern Alpine Ridge at slightly different, yet very close locations near the village of Baceno and the mountain peak Pizzo Tambo. Time series are extracted from the single closest grid points to Baceno (for the 2.2 km simulation) and Pizzo Tambo (for the 12 km simulation) respectively (red points in Fig. 1).

300
To disentangle impacts of the convection parameterization from impacts of spatial averaging, another time series is derived by remapping the 2.2 km simulation conservatively to the 12 km grid and then again extracting the grid point closest to Baceno.
Furthermore, the extracted time series are also aggregated to different time levels in order to differentiate time step dependence in our model from smoothing effects in the precipitation forcing time series.  Zeman et al. (2021). Native and Agg. mean the scale (both temporal and spatial) of the original and the aggregated timeseries respectively.

305
A, numerical convergence is measured with respect to a high-resolution reference simulation (∆z = 1 cm, ∆t = 1s, n z = 200)

Comparison to the Analytical Solution of the Fully Saturated Steady-State in Experiment A
Before the results from the numerical experiments are presented, an analytical solution to the Richards equation is presented for a special case. We consider the case of a fully saturated soil where the volumetric soil moisture corresponds to the total 310 pore volume (θ = η). Furthermore let us assume stationarity (∂θ/∂t = 0), by assuming that the water loss due to ground runoff is always balanced by infiltration from above. Under these assumptions, (4) reduces to Given the boundary condition at the bottom K (z = z b ) = 0, it is straightforward to find the analytical solution for (19). As there is no explicit dependence of K on z, (19) can be solved by separation of variables (∀z ∈ [0, z b ]) and the solution formally 315 reads, Plugging this into the Schlemmer et al. (2018) formulation for ground runoff given by (5), the integral is straightforward to evaluate and yields

320
As Q is given by the Schlemmer et al. (2018) model and K by the soil hydraulic model, the equation is inconsistent, except if the free parameter L −1 g is chosen in an appropriate way which may easily be derived from (19). If one chooses a free drainage boundary condition, the consistency of (19) is naturally given (but only in the absence of a sink, such as e.g. evapotranspiration).
In the case of the (Schlemmer et al., 2018) runoff, it is likely that the flux corrections adjust K to match (21), this is however not included in their analytical setup of the model. The inconsistency arises, because the saturation form of Richards equation 325 is (strictly speaking) not applicable for fully saturated soils (A). Nevertheless, if K is allowed to adjust, ground runoff Q dictates the steady-state solution for K. This is of course simply a consequence of mass conservation -what drips out at the bottom must be refilled at the top in order to sustain the equilibrium. Such a steady-state solution can only be reproduced with a numerical solver, if the chosen scheme is mass-conserving. For any implementation based on flux corrections (here, in particular Schlemmer et al., 2018), (21) implies that the infiltration should be part of the flux correction, because otherwise 330 there will be a discontinuity in K at z = 0.
There is another important consequence of (19): If the solution of K at the surface exceeds K 0 , the equilibrium solution given by (21) cannot be sustained. Such a solution can only exist, if K(z = 0) ≤ K 0 , which implies that we can find a maximum soil depth z b,max for which a stationary solution with a saturated soil is possible, Solutions for for ∆t = 1s with varying grid spacings.
Or -looking at it the other way round -this equation also determines the highest possible position of the water table for a given combination of L −1 g and S ORO . As an example, let L −1 g = 0.4 m −1 , which is the default value for the Schlemmer et al. (2018) model and S ORO = 0.2 which corresponds to a grid point e.g. in a mountainous environment. In this example, we find z b < 12.5 m.

340
The black line in Fig. 4 shows the analytical solution for hydraulic conductivity K for the fully saturated steady-state case.
A time slice at t = 22 h is then extracted from all simulations in experiment AR (Experiment A, regular grids). The numerical solutions found for the different setups reproduce the analytical solution accurately and a selection is shown in Fig. 4. Evidently, the fully saturated steady-state is not of much practical relevance, but the good numerical convergence indicates that the scheme presented in Schlemmer et al. (2018) is indeed mass-conserving irrespective of the chosen grid.

Convergence in Infiltration Experiment A
The convergence during the infiltration process is more relevant for water partitioning at the surface and will therefore be investigated in more detail in this section. Fig. 5 (a) -(c) show time slices of saturation θ f after 1.5, 3 and 5 h from the reference run (∆t = 1 s, ∆z = 1 cm). For all time slices, there is a sharp drop of saturation from θ f = η to θ f = 0.8η, which marks the position of the wetting front. While the position of the wetting front changes with time, its shape remains more 350 or less constant over time, in particular the sharp gradient of saturation is preserved and the propagation velocity is roughly uniform. Apart from the infiltration process that may be seen through the propagation of the wetting from panel (a) to panel (c), there is also an accumulation of water at the bottom of the soil column, which is due to the ground runoff formulation which

Equidistant Grid, AR
The systematic slowdown of the wetting front with increasing time step calls for a more systematic analysis which is presented converge better. Still, the NRMSD values are around 30 % for most model setups with the MVG parameterization, which is considerable. One important reason for the better convergence of the MVG parameterization is that the diffusion close to saturation is much higher than in the RT parameterization. This weakens the sharp gradient around the wetting front and it is therefore easier to resolve it with coarser resolution setups. An experiment with the hydraulic diffusion at saturation 370 D M V G (θ = η) limited to that of the RT parameterization exhibits NRMSD values that are comparable in magnitude to those found for the RT parameterization ( B). Fig. 7 shows NRMSD values for the geometric grid setups. As is the case for regular grids, simulations with the MVG parameterization tend to converge better. Reasonable convergence in the sense that the NRMSD does not exceed 20 % is only achieved 375 for the simulations with b = 0, n z = 200 (which corresponds to the regular 1 cm grid) and for b = 1.01 · 10 −2 , n z = 110 and again only for time steps less or equal to 40 s.

A Closer Look at the Wetting Front
The relatively high NRMSD for long time steps and low resolution grids in Figs. 6 and 7 raises the question whether one can identify an underlying reason for the rather poor numerical convergence. It is illustrative to reconsider Fig. 2 and to again ap-380 preciate the sharp gradient between dry and wet zones in the sand. The formation of such gradients is an inherent characteristic of the dynamics described by Richards equation and it is critical to resolve them as well as possible. The propagation of the wetting front is not solely driven by the gravitational part of the flux K, the diffusive flux D plays an equally important role as shall be explained in the following. The contribution of the two fluxes to the moistening of the soil is shown in Fig. 8 in terms of tendencies, which are given by First, consider the tendencies for the reference run (top row). As is the case for the saturation gradient, the shape (but not the vertical location) of the tendencies is very similar after 1.5 hours and after 6 hours. The diffusive flux acts downstream of the conductive flux, effectively smearing out the saturation gradient slightly. While the contribution of the diffusive flux to 390 the total tendency is smaller than the contribution of the conductive flux, it plays a major role in 'pre-wetting' the soil matrix in front of the tendency maximum. This increase in saturation θ f in advance of the wetting front subsequently enhances K, which is very sensitive to θ f in both (MVG and RT) soil hydraulic models. The anatomy of the wetting front is important in order to understand, why coarse resolution setups do not only fail to capture the sharpness of it but subsequently also fail to reproduce the correct propagation velocity. The thicker the layers underneath the wetting front are, the slower θ f increases.

395
This is illustrated by the broader saturation gradient in a coarse resolution setup (∆z = 10 cm) (Fig. 8, bottom row, blue line).
As there is a strong, nonlinear dependency of D and K on θ f , slower moistening of a layer leads to a delayed acceleration of the fluxes and subsequently to a slowdown of the wetting front with respect to the high resolution setup. This becomes very clear when considering setups with geometric grids as shown in the middle row of Fig. 8. After 1.5 hours, the position and shape of the wetting front is comparable to the wetting front of the reference run, because the layer thicknesses are similar in 400 the uppermost part of the soil. After 6 hours, the shape of the wetting front looks more similar to the one found in the low resolution setup and the position of the wetting front is delayed with respect to the reference run, but not as much as is the case for the coarse resolution setup -indicating a decrease in propagation velocity with increasing layer thickness. The consequence on water partitioning is immediately clear from the position of the wetting front after 6 hours: It is located more than 25 % further up in the low resolution simulation, which of course translates directly to more than 25 % less water in the soil at that 405 moment and it will become more, the longer the infiltration process lasts. While the situation looks better for the telescope setup, it is evident from the considerations in this section that the convergence will worsen over time.    Table 4).
The time series are extracted from individual grid points of the 2.2 km simulation (blue), the 12 km simulation (orange), and the 2.2 km aggregated to the 12 km grid (red). The location of the grid points is indicated in Fig. 1 .

Convergence in Infiltration Experiment B
The rather poor convergence in experiment A is of some concern, but due to the constant high intensity precipitation forcing, the implications might be less severe in a more realistic setup. Time series are extracted from the CRM and CPM simula-410 tions described in section 3.3. An overview can be found in Table 4 and a selection of three time series (B1, B3 and T1) is shown in Fig. 9. The highest intensity precipitation peak is found in the CPM simulation shortly after 00 UTC on May 29, but subsequently the CRM simulation exhibits higher precipitation intensities. The intensities in the CRM simulation aggregated to the 12 km grid are lower than the intensities found for the native grid but the temporal evolution is very similar. This is an expected effect of the spatial averaging. An interesting sidenote concerns the well-known 'flickering' of the convection 415 scheme (see e.g. Guichard et al., 2004) which is visible in the time series of the CPM simulations between 9:00 and 15:00 UTC.
The infiltration process as a response to the timeseries from the convection resolved model for different soil model setups is shown in Fig. 10. For simplicity, the analysis of experiment B will be focused on just four different grid setups with (n z = 200, 110, 32, 8). At first glance, the process seems to be captured in a similar way by all simulations except for the low 420 resolution simulations with n z = 8. In this case, the resolution is simply too coarse to resolve the propagation of the wetting front. For the intermediate resolution case (n z = 32), the propagation of the wetting front is captured, but the accumulation of water in the deep soil is very different from the reference (n z = 200) and high resolution setups (n z = 110), leading to a diverging behaviour of the Schlemmer et al. (2018) model (not shown). This is beyond the scope of this paper as it is model specific and the focus is on the infiltration process which is more similar across different LSMs. At closer look, it can also be seen that the infiltration of water reaches deeper for the ∆t = 20 s than it does for the ∆t = 90 s, while the effect of the input time step ∆t r seems to be negligible.
Still, the time step and grid dependency shown in Fig. 10 is somewhat inconclusive and the consequences on water partitioning are again investigated by comparing surface runoff of the individual simulations. Fig. 11 summarizes the surface runoff 430 response of different simulations to the precipitation forcing given by the the time series in Table 4. Here, accumulated runoff is shown as it is the quantity that will eventually determine soil water availability after rainfall events. The focus is thereby on the first 16 hours as the precipitation amount seems negligible afterwards (Fig. 9). The most striking feature of Fig. 11 are the almost exactly identical surface runoff curves produced by grid setups with n z = 200, n z = 110 and n z = 32, which is why one can identify virtually only two different lines on all panels. The accumulated surface runoff for the low resolution case 435 n z = 8 differs significantly from simulations with other grid setups. It is very important to note that the differences in absolute amounts of water are most pronounced for short time steps and precipitation forcing from high resolution CRM simulations (top row) -the situation found in convection resolving NWP and climate models, where grid setups with eight layers covering roughly two meters of soil (or less) are common (e.g. Heise et al., 2003;Balsamo et al., 2009;Niu et al., 2011;Ducharne, 2016;Reick et al., 2021). The differences between low and high resolution setups are still considerable for simulations with forcing 440 time series from the aggregated CRM (B3 and B4, middle row), but small for simulations with forcing from CPM models (T1, bottom row). As is the case in Fig. 10, the temporal aggregation of rainfall to ∆t s = 90 s has a negligible effect on the simulation's response (first two panels of the top and the middle row). In stark contrast, if the model time step is raised to 90 s, the differences are clearly visible and the effect is most pronounced for the 2.2 km precipitation forcing (compare top row, 2nd and 3rd panel). The systematic increase of surface runoff with increasing time step is evident for all considered forcing time 445 series. Interestingly, the importance of the vertical resolution seems to decrease with increasing time step, indicating that the time step is the dominant factor for convergence.
Under CRM forcing, low resolution simulations with ∆t = 180 s and n z = 8 produce more than twice the amount of surface runoff than high resolution simulations with ∆t = 20 s, n z = 200. If only the effects of vertical resolutions are considered, the differences still amount to more than 50 % of the surface runoff formed in the high resolution simulation. These are 450 considerable amounts of water which are 'lost' in low resolution setups, simply due to the poor numerical convergence.

Conclusions
In this work, we investigated the sensitivity of a common numerical implementation of the Richards Equation with respect to spatial and temporal resolution, and addressed the implications on infiltration and surface runoff processes in land-surface models. In particular, we assessed the effects of the precipitation intensity, time step and vertical discretization on the hydro-   (2018). The vertical discretization was improved by introducing a coordinate transform which allows for a 2nd-order accurate representation of the solution. The model was then applied to investigate the infiltration process. In a first experiment, the model was exposed to constant precipitation rates of 20 mm per hour to investigate the convergence of the infiltration process. 460 We then carried out a second experiment with precipitation time series obtained from convection resolving and convection parameterized simulations to asses the impact of numerical convergence on water partitioning at the surface. The main findings of these experiments are: -Poor convergence for regular grids with a vertical resolution coarser than 4 cm.
-Poor convergence for time steps longer than 40 s.

465
-Better convergence for the Mualem -van Genuchten model than for the Rijtema model.
-Similar situation for unequally-spaced geometric grids, poor convergence for grids with 49 or less vertical layers on a total depth of 2 m.
-Systematic decrease of the propagation velocity of the wetting front with decreasing vertical resolution and increasing time step.

470
-When forced with realistic precipitation time series, the amount of generated surface runoff is more than 50 % higher for low resolution setups than for high resolution setups.
-If long time steps are chosen, the situation is even worse with low resolution setups generating up to 100 % more surface runoff than their high-resolution counterparts.
It is assumed that the strong gradients around the wetting front must be resolved in order to capture the propagation of the 475 wetting front accurately. If the vertical resolution or the time stepping fail to resolve the propagating wetting front, the poor convergence of the infiltration process translates directly to a poor convergence of water partitioning at the surface.
Seven or eight layers covering a two meter soil column is not uncommon in many state of the art LSMs. At these resolutions it might be not appropriate to think of the scheme as an approximation to Richards equation. Instead one may think of it as moving water between boxes of fixed size. This is important as in the latter case, the choice of box size matters a lot for the 480 dynamics of the system and this should be made transparent for model developers and users alike to avoid unwanted effects when the layering is changed.
The detailed analysis illustrates that there are compensating effects at play: while rainfall intensities increase in CRMs, favoring more surface runoff, the employed time step decreases, enabling a better numerical representation of the infiltration process and thus reducing unrealistic runoff.

485
In the full comparison between CPM and CRM simulations, other factors like differing soil types and interception reservoir sizes interfere with the uncertainties arising due to inaccurate numerical approximations of the Richards equation. If the aim is -and it must be -to have similar representations of the terrestrial water cycle across horizontal resolutions, model tuning is In order to appreciate the issue related to the saturation form Richards equation, we revisit its derivation from (3) to (4), where the hydraulic diffusivity is introduced as D(θ) = K(θ)∂h/∂θ. It must be stressed that the introduction of D is not motivated by the physics of the system, but solely in order to bring the equation to a form for which solution methods readily existed at the time (namely the diffusion equation). In soil physics, the expression C(θ) = ∂θ/∂h is referred to as specific water capacity (see Hillel, 2012, p.158) and thus hydraulic diffusivity is the ratio of hydraulic conductivity to specific water capacity. Clearly, 505 in the saturated part of the soil, specific water capacity must vanish and as K(θ) is finite, hydraulic diffusivity diverges, Therefore the chain rule expansion ∂h/∂z = ∂h/∂θ · ∂θ/∂z is mathematically invalid for fully saturated soils as it is an expansion with "∞ · 0". This implies that one should not use the saturation form of Richards equation when the soil is expected to reach saturation. Formally, this is shown by Gilding (1991) who derives integrability conditions for D under which wetting 510 fronts with finite propagation speed may form, at least one of which is not met by (A1). Therefore the saturation form is neither suitable to treat the infiltration problem, nor is its use adequate in the presence of water tables. This is also the underlying reason for the necessity of limiting hydraulic diffusivity in the MVG parameterization.
Apart from the issue with diverging hydraulic diffusivity pointed out in the beginning of this section, there is another closely related problematic aspect in the saturation form: Close to the water table, discontinuities in θ may appear and thus render θ 515 not differentiable in space. In contrast, the suction head h is continuous and differentiable even in the transition region from vadose zone to water table (Farthing and Ogden, 2017).
Despite the issues explained in this appendix, it should be mentioned that the saturation form is used successfully in LSMs, even if its use is mathematically not fully appropriate. In general, the product D(θ) · ∂θ/∂z converges (otherwise the mixed to finite values. The alternative is to either use the mixed form (3) or a fully head based form which can be derived from (3) by using the chain rule to expand ∂θ/∂t = ∂θ/∂h · ∂h/∂t and plugging in the definition of the specific water capacity to obtain The fully head-based form is usually favored by soil physicists. The reason is that if one accounts for the compressibility of 525 water, the left hand side of (3) would contain an additional term in which the head is appearing. The numerical solution of the mixed form is thus more complicated as one needs to advance both θ and h in time (Tocci et al., 1997). In the absence of compressibility effects, the mixed form is however appropriate and solution methods exist in the soil physics literature (De Rooij, 2010).