Articles | Volume 27, issue 13
Research article
04 Jul 2023
Research article |  | 04 Jul 2023

Assessment of the interactions between soil–biosphere–atmosphere (ISBA) land surface model soil hydrology, using four closed-form soil water relationships and several lysimeters

Antoine Sobaga, Bertrand Decharme, Florence Habets, Christine Delire, Noële Enjelvin, Paul-Olivier Redon, Pierre Faure-Catteloin, and Patrick Le Moigne

Soil water drainage is the main source of groundwater recharge and river flow. It is therefore a key process for water resource management. In this study, we evaluate the soil hydrology and the soil water drainage, simulated by the interactions between soil–biosphere–atmosphere (ISBA) land surface model currently used for hydrological applications from the watershed scale to the global scale, where parameters are generally not calibrated. This evaluation is done using seven lysimeters from two long-term model approach sites measuring hourly water dynamics between 2009 and 2019 in northeastern France. These 2 m depth lysimeters are filled with different soil types and are either maintained as bare soil or covered with vegetation. Four closed-form equations describing soil water retention and hydraulic conductivity functions are tested, namely the commonly used equations from Brooks and Corey (1966) and van Genuchten (1980), a combination of the van Genuchten (1980) soil water retention function with the Brooks and Corey (1966) unsaturated hydraulic conductivity function, and, for the very first time in a land surface model (LSM), a modified version of the van Genuchten (1980) equations, with a new hydraulic conductivity curve proposed by Iden et al. (2015). The results indicate good performance by ISBA with the different closure equations in terms of soil volumetric water content and water mass. The drained flow at the bottom of the lysimeter is well simulated, using Brooks and Corey (1966), while some weaknesses appear with van Genuchten (1980) due to the abrupt shape near the saturation of its hydraulic conductivity function. The mixed form or the new van Genuchten (1980) hydraulic conductivity function from Iden et al. (2015) allows the solving of this problem and even improves the simulation of the drainage dynamic, especially for intense drainage events. The study also highlights the importance of the vertical heterogeneity of the soil hydrodynamic parameters to correctly simulate the drainage dynamic, in addition to the primary influence of the parameters characterizing the shape of the soil water retention function.

1 Introduction

Drainage water is the portion of precipitation that flows through the first metres of soil. As it has escaped evapotranspiration, it is the main diffuse source of groundwater recharge (Moeck et al.2020) and is of crucial importance for estimating the evolution of an aquifer (Bredehoeft2002; Döll and Fiedler2008). Even where there is no aquifer, this water flux can contribute to the baseflow. Despite their importance, direct observations of drainage water are still rare compared to direct observations of river flow or groundwater level.

Indirect methods based on an analysis of baseflow (Meyboom1961) or a variation in the piezometric level, for instance, using a water table fluctuation method (Healy and Cook2002), can be applied to estimate groundwater recharge. However, these methods cannot separate the different components of the recharge, such as the exchange between surface water and groundwater (Brunner et al.2017; Keshavarzi et al.2016), exchanges between several aquifer layers basins (Tavakoly et al.2019), and drainage water. Direct measurement at the local scale and at a high frequency can be made in situ, using lysimeter columns, which isolate a small volume of soil from the lateral flow and collect drainage water. Most of the time, the lysimeters are disconnected from the soil and avoid the capillary rise from the deeper soil that can have an important influence locally (Vergnes et al.2014; Maxwell and Condon2016). From such observations, drainage water is known to have large variations in space, due to changes in soil texture and structure (Vereecken et al.2019; Moeck et al.2020).

Due to the limited number of observations and the difficulty around indirectly quantifying it, the estimation of the groundwater recharge is 1 of the 23 unsolved problems in hydrology (Blöschl et al.2019). The simulation of recharge in hydrological models can vary significantly, from simplified reservoir approaches to physically based models. The widely used reservoir approach (Alcamo et al.2003; Harbaugh2005) has the advantage of limiting the number of parameters to be calibrated and reducing the numerical cost of simulations. Physically based approaches are more complex and are commonly used in land surface models (LSMs) in which vertical water and energy balances between the land surface and the atmosphere can be calculated (De Rosnay et al.2003; Blyth et al.2010; Boone et al.2000).

Today, LSMs are widely used in hydrological applications in order to study the regional and global water cycle, predict streamflow, and inform water resource management (Haddeland et al.2011; Vincendon et al.2016; Schellekens et al.2017; Gelati et al.2018; Le Moigne et al.2020; Vergnes et al.2020; Muñoz-Sabater et al.2021; Rummler et al.2022). In this context, the challenge of LSMs is finding a compromise between a simple application and an application that is powerful enough to reproduce the full water cycle. For example, in the unsaturated zone, hydrodynamic parameters are generally uncalibrated and are estimated with soil properties (Decharme et al.2011, 2019; Le Moigne et al.2020). However, LSMs may have difficulty estimating the dynamics of groundwater recharge, particularly during intense precipitation events. Vereecken et al. (2019) suggested a number of directions for improvement, namely introducing more physical processes such as preferential flow along the roots and macropores, improving the representation of the soil and vegetation and of soil parameters, and improving the spatiotemporal distribution of precipitation. In LSMs, the Richards (1931) equation is used to describe the flow of water in the porous soil due to the actions of capillarity and gravity. Mainly used in the hydrology community, this equation has been widely criticized, in particular for the overestimation of the effect of capillarity (Nimmo2010; Farthing and Ogden2017).

In the Richards (1931) equation, two closed-form equations, namely Brooks and Corey (1966) and van Genuchten (1980), hereafter BC66 and VG80, respectively, are often used to represent the variations in the volumetric water content with the matrix potential and the hydraulic conductivity in the unsaturated zone. Historically, the closed-form equations from BC66 are mostly used by the atmospheric community in theirs LSMs, while the ones from VG80 are mainly used by hydrologists. BC66 proposed the simple analytical power curves of soil water retention and hydraulic conductivity, based on North American soil observations. However, they do not include an inflection point close to saturation and are thus not derivable, leading to problems at the connection with the saturated zone. The consequence is a deviation near the saturation of volumetric water content. VG80 proposed an improvement of the BC66 soil water retention curve close to water saturation, using more complex analytical forms that are based on European soil observations. However, as summarized by Iden et al. (2015), when such an S-shaped function is used to described fine-textured soil or a soil with a wide pore size distribution, the VG80 hydraulic conductivity curve exhibits an abrupt drop at the transition from saturated conditions to unsaturated conditions. This may lead to underestimations of the hydraulic conductivity for very wet conditions, i.e. near saturation, and impact the stability and the accuracy of the numerical solver (Vogel et al.2000; Ippisch et al.2006). The nonlinear form also increases the cost of the numerical solvers (van Genuchten1980; Vogel et al.2000; Ippisch et al.2006; Dourado Neto et al.2011).

To smooth the artefacts associated with the VG80 hydraulic conductivity curve, some studies suggest different solutions, such as shifting the entire pore size distribution with an air entry value (Kosugi1994), truncating the pore size distribution by introducing a maximum pore radius on the VG80 hydraulic conductivity curve (Schaap and Van Genuchten2006; Iden et al.2015), or putting an explicit air entry pressure in the VG80 model (Fuentes et al.1992; Braud et al.1995; Vogel et al.2000; Ippisch et al.2006). Fuentes et al. (1992) or Braud et al. (1995) proposed combining the VG80 soil water retention function with the BC66 hydraulic conductivity curve. Such a combination should improve the simulation when the soil water is close to saturation, while preserving the simplicity and numerical stability of the BC66 relationships. It benefits from the many methods that link soil parameters from one relationship to another (van Genuchten1980; Lenhard et al.1989). Recently, an elegant parametrization was proposed by Iden et al. (2015), which allows us to keep the VG80 hydraulic conductivity curve by introducing the concept of maximum pore radius, as suggested by Durner (1994).

The goal of this study is to use a state-of-the-art LSM to assess the benefits of the BC66 and VG80 relationship and two of its derivatives in reproducing soil water mass, volumetric water content, and drainage flux observed in seven lysimeters during more than 5 years. The two alternative approaches are (i) the one proposed by Fuentes et al. (1992) and Braud et al. (1995) that combines the VG80 soil water retention curve with the BC66 unsaturated hydraulic conductivity curve and (ii) the new parameterization of the VG80 hydraulic conductivity curve from Iden et al. (2015) that has never been used in an LSM, to the best of our knowledge. The soil hydrodynamic parameters are directly derived from observations and compared with several pedotransfer functions (hereafter PTFs) commonly used by LSMs. The LSM used is the multi-layer diffusion version of the interactions between soil–biosphere–atmosphere (ISBA) that directly solves the mixed form of the Richards (1931) equation (Boone et al.2000; Decharme et al.2011). The experimental protocol, including descriptions of both the lysimeter data and the ISBA model, is given in Sect. 2. Soil parameter estimation based on lysimeter observations is described in Sect. 3, while the main results of each model approach are presented in Sect. 4. Finally, a general discussion and the main conclusions are given Sect. 5.

2 Experimental protocol

2.1 Data

Seven lysimeters from two French experimental sites located in northeastern France are used in this study. There are three lysimeters from the GISFI experimental station (French Scientific Interest Group – Industrial Wasteland) at Homécourt (4921 N, 599 E; altitude 430 m), with a data record from 2009 to 2016 (hereafter G1, G2, G3, and G4), and three lysimeters from OPE experimental station (Perennial Environmental Observatory) close to Osne-le-Val (485092 N, 52119 E; altitude 224 m), with a data record from 2014 to 2019 (hereafter O1, O2, and O3). These two sites are separated by a distance of 97 km. Their soils are classified according to the World Reference Base for Soil Resources (WRB; Hazelton and Murphy2016). For the lysimeters, the rain intercepted by the vegetation is assumed to evaporate or to eventually reach the ground, as this is the case for all natural surfaces with a vegetation canopy.

The GISFI experimental station focuses on the understanding of pollution evolution and the development of decontamination solutions (Lemaire et al.2019; Huot et al.2015; Rees et al.2020; Ouvrard et al.2011). It participates in the collaborative study of TEMPOL (Observation sur le long TErMe de sols POLlués) with the German observatory infrastructure of TERENO (Terrestrial Environmental Observatories; Zacharias et al.2011) in order to study in situ the transfer of pollutants (Leyval et al.2018). Three of the GISFI lysimeters (G1, G2, and G3) contain rebuilt soil of the Spolic toxic Technosol, which was sampled in Neuves-Maisons, an industrial wasteland of a coking plant with contamination (polycyclic aromatic hydrocarbons (PAHs) and hydrocarbons; Monserie et al.2009; Ouvrard et al.2011; Sterckeman et al.2000). These three lysimeters were filled in September 2007 and the soil material was gradually and manually compacted every 0.3 m to reach a given bulk density. Lysimeters G1 and G2 were maintained as bare soil, while lysimeter G3 was covered by vegetation (alfalfa). Lysimeter G4 was filled by a monolith of Cambisol from Noyelles-Godault and covered by grass (Table 1).

Table 1Description of each lysimeter, including the filling method, soil type, vegetation cover, number of texture observations, and textures (in percent of clay, sand, and silt) at different depths.

* Hypereutric Cambisol.

Download Print Version | Download XLSX

The main objective of the OPE site is to describe the environment before and after the construction of the surface facilities of a deep radioactive waste repository and to follow its evolution. This site is part of the OZCAR (Observatoires de la Zone Critique: Application et Recherche) French research infrastructure dedicated to the observation and study of the critical zone (Gaillardet et al.2018; Braud et al.2020). The lysimeters were filled by monoliths taken close to the OPE site; lysimeters O1 and O3 were filled with hypereutric Cambisol, with different layers of limestone (more or less cracked), and lysimeter O2 contains a Cambisol. These three lysimeters had sparse vegetation composed of C3-type grass (Table 1). Monolith filled lysimeters preserve original soil horizons and are therefore better representative of watershed conditions.

All of the lysimeters are weighable cylinders, with a depth of 2 m and an area of 1 m2. No suction is imposed at the bottom of each lysimeter, contrary to the ones used for TERENO-SOILCan (Putz et al.2016). They are equipped with suction and temperature probes, in addition to time domain reflectometry (TDR) probes to measure the water content, at 0.50, 1.0, and 1.5 m depth, with an additional measurement at 0.2 m for the OPE lysimeters. The weight is measured continuously with a 0.1 kg precision, giving an idea of the time variations in the total soil water volume in the column, and a tipping counter measures the drainage water at the bottom. Data are continuously monitored on an hourly basis using a data logger. On the GISFI site, TDR probes are TRIME-PICO 32 sensors with internal TDR electronics. They are set horizontally and record the water content (in cm3 cm−3) on an hourly basis. The calibration is performed on two measurements, namely one in a dry condition and one in a water-saturated condition. At the OPE site, the soil moisture sensors used (UMP-1; Umwelt-Geräte-Technik GmbH) are based on a frequency domain reflectometry (FDR) method and measure the local change in dielectric permittivity. The tipping bucket resolution is 0.1 mm h−1 at the two sites.

The sets of atmospheric forcing variables (wind speed, precipitation rate, short-wave incident radiation, air temperature, air humidity, and atmospheric pressure) used to force the ISBA LSM are observed in situ at an hourly time step by two local meteorological stations (one at OPE and one at GISFI). The atmospheric forcing is assumed to be identical for all of the lysimeters at each site. Annual precipitation is 20 % higher at the OPE site (876 mm yr−1) than at the GISFI experimental station (727 mm yr−1). Long-wave radiation is derived from atmospheric conditions, using the equation of Prata (1996). Atmospheric data gaps are filled by regressions on available data, using two neighbouring meteorologic stations. The gaps represent up to 12 % of the observations for the GISFI site.

All soils for all of the lysimeters are described in Table 1, while all available observations and their characteristics are summarized in Table S2.

2.2 ISBA model

The ISBA LSM was originally scheduled for use in numerical weather prediction and climate models. During the last few decades, ISBA has evolved from a simple bucket force–restore model (Noilhan and Planton1989) to a more explicit soil multi-layer diffusion scheme (Boone et al.2000; Decharme et al.2011) coupled with a snow multi-layer scheme of intermediate complexity (Boone and Etchevers2001; Decharme et al.2016). ISBA is currently used in hydrological applications, especially to estimate groundwater recharge when it is associated with hydrologic models at both regional (Le Moigne et al.2020; Tavakoly et al.2019; Vergnes et al.2020) and global scales (Vergnes and Decharme2012; Decharme et al.2019). Several studies showed the good performance of ISBA in simulating field sites (Calvet et al.1998; Boone et al.2000; Decharme et al.2011; Joetzjer et al.2015; Garrigues et al.2015, 2018; Morel et al.2019).

The surface temperature evolves with the heat storage in the soil–vegetation composite and the thermal gradient between surface and the other layers (Boone et al.2000; Decharme et al.2011). At depth, the heat transfer is described by the use of the one-dimensional Fourier law. The soil heat capacity is the sum of the water heat capacity and the heat capacity of the soil. The soil thermal conductivity is a function of volumetric water content and porosity. Freezing due to water phase changes in the soil can also be computed by taking into account the effect of the sublimation of frost and the isolation of vegetation at the surface (Decharme et al.2016).

ISBA simulates the land surface evapotranspiration as being the sum of the bare-soil evaporation, the soil freezing sublimation, the plant transpiration, the direct evaporation of the precipitation intercepted by the plant canopy, and the snow sublimation (Noilhan and Planton1989). Water for bare-soil evaporation is drawn from the first layer of the soil. This soil evaporation is weighted by the relative humidity of this superficial layer (Mahfouf and Noilhan1991). This relative humidity evolved nonlinearly with the superficial water content, potentially allowing the moisture content of the soil evaporation to be greater than the usual water content at the field capacity, which is specified as the matric potential at −0.33 bar. Plant transpiration is proportional to a stomatal conductance computed using an interactive vegetation scheme. The water used for transpiration is removed throughout the root zone in which the roots are asymptotically distributed, according to Jackson et al. (1996). Plant transpiration stops when the root zone water content is below the usual water content at wilting point, corresponding to a matric potential of −15 bar.

The interactive vegetation scheme of ISBA activated in this study represents 16 broad types of vegetation. The scheme represents plant photosynthesis and respiration, plant growth, and decay. The simulated stomatal conductance depends on photosynthesis and controls transpiration. The vegetation canopy is characterized by the leaf area index (LAI), which results directly from the carbon balance computation of the leaves. The simulated LAI varies during the year, according to photosynthetic activity, respiration, and decay. In spring, for instance, when photosynthetic activity is high due to high solar radiation and sufficient volumetric water content, the LAI increases. LAI affects the transpiration and the evaporation of water on the canopy associated with intercepted rain or dew deposition. Photosynthesis and respiration are parameterized, according to the semi-empirical model of Goudriaan et al. (1985), and implemented by Jacobs et al. (1996) and Calvet et al. (1998). Plant growth and decay are based on Calvet and Soussana (2001), Gibelin et al. (2006), and Joetzjer et al. (2015). A complete description of the carbon cycle in vegetation can be found in Delire et al. (2020).

In ISBA, the intercepted rain by vegetation is fully represented and based on a simple rainfall interception scheme (Noilhan and Planton1989; Noilhan and Mahfouf1996). The interception reservoir is fed by the intercepted rain by vegetation. When this reservoir is larger than its maximum value (equal to 0.2xLAI), then the dripping from the vegetation is computed as a simple mass balance. The direct evaporation of the vegetation is drawn from this interception reservoir and depends on the fraction of the foliage covered by intercepted water, as proposed by Deardorff (1978).

The ISBA soil hydrology uses the following mixed form of the Richards (1931) equation to describe the water mass transfer within the soil via Darcy's law:

(1) ω ( z ) t = z - k ( z ) ψ ( z ) z + 1 + S ( z ) ρ ω ,

where ω (m3 m−3) is the volumetric water content at each depth z (m), ψ (m) ia the water pressure head, and k (m s−1) is the soil hydraulic conductivity. S (kg m−3 s−1) is the soil water source/sink term related to water withdrawn by evapotranspiration in each layer. A complete description of the model equations used to simulate water transport can be found in Boone et al. (2000) and Decharme et al. (2011, 2019). This soil hydrology is solved numerically, using a Crank–Nicolson time scheme, where the flux term is linearized via a first‐order Taylor series expansion. The resulting linear set of diffusion equations can be cast in a tridiagonal form and solved quickly (see Sect. S1 in the Supplement for more details). The soil discretization is adapted to the lysimeter depth and to the measurement depths for which 14 soil layers are used at 0.01, 0.04, 0.1, 0.2, 0.4, 0.6, 1.0, 1.2, 1.4, 1.6, 1.8, 1.95, and 2 m depths. ISBA uses the soil infiltration at the surface as upper boundary condition, neglecting the surface runoff when simulating local field or lysimeter sites. This soil infiltration is the flux of water reaching the soil surface (i.e. the sum of the precipitation not intercepted by the canopy, the dripping from the vegetation, and the snowmelt). At the bottom of the soil column, a fine layer of 5 cm is used to simulate a seepage face lower boundary condition (LBC) instead of the usual free-drainage LBC used in ISBA (see Sect. S1.3). Such a seepage face LBC is recommended to simulate the lysimeters (Séré et al.2012; Tifafi et al.2017). Finally, to compute the initial conditions, a spin-up of 1 year is done, using the first year of each forcing dataset in order to ensure an adequate numerical equilibrium for soil water.

2.3 Model approaches

In this study, we evaluate four closed‐form equations that link the volumetric water content, soil matric potential, and hydraulic conductivity, assuming a residual water content (ωr) equal to 0. All simulations use the same ISBA configuration and are denoted as follows:

  • BC66. This first model approach uses the standard version of ISBA, i.e. the closed‐form equations of Brooks and Corey (1966), where the soil water retention function, ψ(ω), and the hydraulic conductivity function, k(ψ), are given by the following:

    (2) ψ ( ω ) = ψ sat ω ω sat - b k ( ψ ) = k sat ψ ψ sat - 2 b + 3 b ,

    where b (–) represents the dimensionless shape parameter of the soil water retention function, while ψsat (m) and ksat (m s−1) are the soil matrix potential and hydraulic conductivity at saturation, respectively.

  • VG80. This second model approach uses the following closed‐form equations from van Genuchten (1980), using the Mualem (1976) theory:

    (3) ψ ( ω ) = - 1 α ω ω sat - 1 / m - 1 1 / n with m = 1 - 1 n k ( ψ ) = k sat S l [ 1 - ( 1 - S 1 / m ) m ] 2 with S = 1 + | α ψ | n - m ,

    where α (m−1) is the inflection point at which the slope of the soil water retention function (dω/dψ) reaches its maximum value, n (–) is a dimensionless coefficient that characterizes the shape of the retention curve, and l is the Mualem (1976) dimensionless parameter that determines the slope of the hydraulic conductivity function.

  • BCVG. This third model approach (see Fuentes et al.1992, or Braud et al.1995) uses a combination of the ψ(ω) function from VG80 with the k(ψ) function from BC66, using the Burdine (1953) theory as follows:

    (4) ψ ( ω ) = - 1 α ω ω sat - 1 / m b - 1 1 / n b with m b = 1 - 2 n b k ( ψ ) = k sat S b 2 λ b + 3 with λ b = n b - 2 - 1 and S b = 1 + | α ψ | n b - m b ,

    where nb (–) is the dimensionless coefficient that characterizes the shape of the retention curve using Burdine (1953) theory.

  • VGc. This last approach uses the usual ψ(ω) function from VG80 in Eq. (3) and the corrected form of the k(ψ) function from Iden et al. (2015), using the Mualem (1976) theory as follows:

    (5) k ( ψ ) = k sat S l Γ 2 1 - 1 - S 1 / m m 2 ψ ψ c 1 - 1 - S c 1 / m m + S - S c | α ψ c | 2 ψ > ψ c ,

    with Sc=[1+|αψc|n]-m and Γ=1-(1-Sc1/m)m+1-Sc|αψc|, where ψc (m) is the value of the maximum suction near saturation. It is fixed at −0.01 m, as suggested by Iden et al. (2015), which is equivalent to use a maximum pore size only in the capillary bundle model.

Additional sensitivity tests were made (see Sect. 5) first to compare the impact of prescribing homogeneous soil profiles, which are classically derived from the use of pedological soil maps, versus heterogeneous soil profiles, which are generally observed in situ. Second, we test the performances of using soil parameter values derived from the usual PTFs (Clapp and Hornberger1978; Cosby et al.1984; Carsel and Parrish1988; Wösten and van Genuchten1988; Vereecken et al.1989; Weynants et al.2009) instead of using soil parameter values observed in situ. Finally, we also assess the impact of using the seepage face LBC compared to the usual free-drainage LBC used in ISBA over natural sites and in regional- and global-scale applications.

3 Estimation of parameters

3.1 Soil parameters

The rich datasets collected by the lysimeters (hourly resolution with 3220 observations on average by lysimeters and by depth) allow us to estimate the soil hydrodynamic parameters, such as in previous studies (Brooks and Corey1966; van Genuchten1980; Carsel and Parrish1988). For instance, Fig. 1 plots the observed soil water retention curves, namely the volumetric water content, ω, as a function of the logarithm of the absolute value of the soil matrix potential, ψ (dots with different colours for each depth) for lysimeters G2 and O1 (other lysimeters are given in Figs. S1 and S2 in the Supplement). To remove the effect of hysteresis on the observed soil water retention curves (i.e. the difference in soil water retention curve measured under wetting and drying process; Haines1930), we have averaged the ω values for each ψ values. Figure 1 reveals an important heterogeneity with depth. From these observed soil water retention curves, ωsat, ψsat, α, n, nb, and b can be derived at each depth of each soil column (for example, b is the slope of this relation). To this end, we use an objective least squares function, which minimizes the sum of the squares of the deviations. This function allows us to maximize the likelihood with a normal distribution (nls function in RStudio). The soil parameter estimates allow us to plot soil water retention curves from BC66 and VG80 (red and blue lines, respectively) in Fig. 1, showing a better fit for VG80 than BC66.

Figure 1Soil water retention curves after removing the effect of hysteresis, with the volumetric water content (ω) and logarithm of the absolute value of the soil matric potential (ψ) for lysimeters O1 and G2. Observations at 0.2, 0.5, 1.0, and 1.5 m depth are represented by orange, aquamarine, grey, and purple diamonds, respectively, and estimations by red and blue dots for BC66 and VG80, respectively. The dashed lines are the estimated values via the observations for the water content at saturation and matric potential.


Next, we use two steps to evaluate this method and the soil parameter estimates. First, we successfully check that (1) the estimated ωsat is consistent with the 99th percentile of the observed soil volumetric water content at each depth and for each lysimeter and that (2) the estimated n, nb, and b or α and ψsat verify the following well-known relationships defined by van Genuchten (1980), as follows:

(6) n 1 + b - 1 for Mualem theory n b 2 + b - 1 for Burdine theory α | ψ sat | - 1 .

Although several authors proposed more complex relationships that allow the preservation of the capillary effects (Lenhard et al.1989; Morel-Seytoux et al.1996; Leij et al.2005; Sommer and Stöckle2010), a comparison of some of these relationships by Ma et al. (1999) has shown that the van Genuchten (1980) and Lenhard et al. (1989) relationships give better results. It is interesting to note that the simple relationship nbn+1 is usually true for our soils (see Table S2). This justifies our choice to not represent the ω(ψ) relationship from BCVG (Eq. 4) in Fig. 1 (and thus in Figs. S1 and S2) because it is similar to the ω(ψ) relationship from VG80 (Eq. 3). Second, we compare our estimates to an alternative method to derive soil parameters from observed soil water retention curves based on the package SoilHyP (Dettmann et al.2022), which uses the shuffled complex evolution (SCE) optimization. Soil parameter values found with this method are very close to our estimates, as shown by the low mean relative bias between both estimates (5 % for b, 1 % for n, 2 % for ωsat, and 4 % for ψsat). Note that, for α, even if the mean relative bias seems to be not negligible (11 %), both estimates still remain very close (20.64 m−1 on average for SoilHyP and 18.46 m−1 for this study).

Estimations of the soil hydraulic conductivity function parameters (ksat and l) are more difficult because of the lack of observed hydraulic conductivity profiles in each lysimeter. When the deepest part of each lysimeter is close to saturation, we assume that the soil hydraulic conductivity at saturation, ksat, can be considered equal to the observed hourly drainage at 2 m depth. ksat is thus derived empirically from the 99.9 percentile of the hourly drainage distributions. No assumption is made on a hypothetical ksat profile with depth; i.e. ksat is taken to be homogeneous with depth in each lysimeter. For VG80, a particular attention is also given to the very low n values for which Eq. (3) is numerically unstable. In this study, we found that our simulations are numerically stable when n≥1.1, so the limit of n is fixed at n=1.1. As already said, BCVG (Eq. 4) and VGc (Eq. 5) allow us to cancel this limitation. The l parameter in Eq. (3) for VG80 and Eq. (5) for VGc is estimated with a simple calibration via the ISBA sensitivity tests, with l ranging from −5 to 5. Better scores are obtained for l equal to 0.5, a classical value, for the OPE experimental station. For the GISFI experimental station, better scores are obtained for l equal to 0.5 (G4), −2 (G1 and G2), and −5 (G3), which remains consistent with the literature (Wösten and van Genuchten1988; Wösten et al.2001; Schaap et al.2001). We choose l to be constant over the entire soil profile for each lysimeter.

Figure 2 gives an example for lysimeters O1 and G2 of the relative soil hydraulic conductivity, k/ksat, as a function of the soil water content actual saturation, ω/ωsat, derived from the four closed-form equations described previously (other lysimeters are given in Fig. S3). The VG80 hydraulic conductivity exhibits an abrupt drop at the transition from saturated conditions to unsaturated conditions, contrary to BC66. This behaviour is well known for fine-textured soil as in our lysimeters (Fuentes et al.1992; Vogel et al.2000; Ippisch et al.2006; Iden et al.2015). This explains why the VG80 hydraulic conductivity function with small n values is strongly unstable near saturation. A comparison with the observed hourly drainage water at 2 m depth (reduced to ksat) versus the actual saturation at 1.5 m depth seems to confirm that VG80 certainly underestimates the actual hydraulic conductivity near saturation for fine-textured soil, as underlined by previous studies (Vogel et al.2000; Ippisch et al.2006; Iden et al.2015).

Figure 2Near-saturation estimates of the relative soil hydraulic conductivity, k/ksat, as a function of the soil water content's actual saturation, ω/ωsat, for lysimeters O1 and G2 at 1.5 m depth. BC66 (red), VG80 (blue), BCVG (green), and VGc (orange) are computed using parameter estimates (Table S2) into Eqs. (2)–(5), respectively. The dots represent the observed hourly drainage water at 2 m depth (reduced to ksat) versus the actual saturation at 1.5 m depth.


All parameter estimates are presented in Fig. 3 and in Table S2. Parameters vary greatly with depth (especially b, n, ψsat, and α) and, to a lesser extent, between lysimeters. nb is not shown in Fig. 3 because it varies as n does. Although these heterogeneities can be observed on field sites, they could also be increased by factors like the compaction and structuration generally present in lysimeters (Weihermüller et al.2007; Séré et al.2012). The largest differences are found for b (and thus n). For all of the lysimeters, except G1 and G4, b increases with depth (and, conversely, n decreases). For the OPE lysimeters, b starts with 8 at 0.2 m depth, doubles twice (first at 0.5 m and then at 1 m depth), and then reaches 50 at 1.5 m. These variations are less pronounced for the GISFI lysimeters, with only a variation with depth by a factor of 2 for G2. The absolute value of ψsat decreases with depth (and, conversely, α increases), especially for the OPE lysimeters. Lysimeters were filled with preserved soil columns at the OPE and manually at the GISFI. Repacked soil columns are generally recognized as being more homogenous than undisturbed monoliths (Weihermüller et al.2007; Carrick et al.2017). The fact that the GISFI lysimeters are filled in a manual way can thus explain the weaker soil heterogeneities with depth, although they have been filled to preserve the bulk density of the soil.

Figure 3From left to right and top to bottom: the hydraulic conductivity at saturation, ksat (106 m s−1), volumetric water content at saturation ωsat (m3 m−3), b and n, matric potential at saturation, ψsat (m), and α (m). Estimations from in situ measurements are represented by triangles at 0.2, 0.5, 1, and 1.5 m (orange, aquamarine, grey, and purple, respectively), and their mean (homogeneous) values are represented by crosses (×). The values derived from the six pedotransfer function are shown by a box plot presenting the median, 25 %, and 75 % quantiles.


When comparing our parameter estimates to the PTFs usually used in LSMs (Clapp and Hornberger1978; Cosby et al.1984; Carsel and Parrish1988; Wösten and van Genuchten1988; Vereecken et al.1989; Weynants et al.2009), we find that b (and thus n) from the in situ estimates are quite different from those determined by the PTFs (box plots). While b ranges from 5 to 12 with the usual PTFs, our estimates range from 8 to 50. ψsat (and thus α) estimates are slightly different from those determined by usual PTFs. For other parameters, our estimates and the usual PTFs give similar values.

3.2 Vegetation parameters

For the lysimeters covered by vegetation (O1, O2, O3, G3, and G4), two additional soil parameters must be determined. The field capacity ωfc and the wilting point ωwilt are computed via matrix pressure at −0.33 and −15 bar, respectively. The root depths have also been determined. Although the rooting depth was not measured at the sites, it is possible to derive it from the observations of the volumetric water content profile; if the volumetric water content presents a slow decrease in summer at a given depth, then it is considered that the roots have not yet reached this depth. The root depth is thus fixed at 2 m for lysimeters G3 and O2 and varies for lysimeters G4, O1, and O3. From 2009 to 2013, the root depth in lysimeter G4 reached 1.5 m; but after June 2013, seeding and harvesting were carried out every year, limiting roots whose development never reached below 0.4 m depth. In lysimeters O1 and O3, the root depth reached 0.8 m from 2014 to 2018 and 2 m thereafter. Standard ISBA values of maximum photosynthesis, leaf nitrogen content, and specific leaf area were used for the lysimeters covered by grass, but specific values were derived form the TRY database (Kattge et al.2020) for alfalfa, as this crop is not a standard vegetation type in ISBA.

As there are no measurements of LAI for these lysimeters, the simulation of the LAI can only be compared to the literature. As expected, the simulated LAI is minimal between December and April and maximal in August, and there are variabilities between years (Fig. S4). For soils with grass (G4, O1, O2, and O3), the maximum LAI varies between 3 and 5, and the mean annual LAI varies between 1.3 and 2. These values are similar to those found in other studies (Calvet2000; Darvishzadeh et al.2008). LAI for alfalfa cover (G3) is larger (LAImax=7.5; LAImean=4.8), which is expected for such well-developed vegetation (Wolf et al.1976; Wafa et al.2018).

4 Results

Here, we present the main results for each model approach described in Sect. 2.3 in terms of simulated soil water and drainage water dynamics, water budget, and intense drainage water events. We used four skill scores to assess the ability of each model approach to reproduce the observation. These include the overall bias between the simulation and observation, the centred root mean square error (CRMSE) computed by subtracting the simulated and observed annual means from their respective time series before computing a standard root mean square error, the coefficient of determination (R2) or the Pearson correlation (R), and the Nash and Sutcliffe (1970) efficiency criterion (NASH) that determines the relative magnitude of the residual variance compared to the measured data variance. These scores are summarized in Appendix A. Note that the simulated soil temperatures have also been studied and analysed. All model approaches have shown that good skill scores (R2>0.9) underlie the ability of the ISBA LSM to reproduce observed soil temperatures. Because there are no significant differences between the four model approaches compared to the observations, these results are not presented in this study.

4.1 Soil water dynamics

The dataset allows the assessment of the evolution of the total soil water mass derived from the mass of the lysimeter, of the soil volumetric water content at several depths, and of the drainage water flux at the bottom of the lysimeters. In the following analysis, periods during which the meteorological forcing is reconstructed or when data are of poor quality (Table S1 in the Supplement) are not taken into account in the score computation. These periods are short, except for the drainage water of lysimeter G3 (23 % of the duration).

4.1.1 Water mass time variations

As there are no observations of the weight of the dry soil in each lysimeter that can serve as a reference, we evaluate the ability of ISBA to simulate the temporal variation in the water mass around the mean and not around the absolute mass of each lysimeter. These variations are assumed to be equal to the time variations in the total mass of the lysimeter, neglecting the seasonal variations in the vegetation mass. Time variations in the water masses are presented in Fig. 4 for the seven lysimeters at 1 h time step. All model approaches (BC66, VG80, BCVG, and VGc) are represented.

Figure 4Hourly time series of the total water mass variations (kg) around the observed or simulated means from GISFI (G1, G2, G3, and G4) and OPE (O1, O2, and O3) lysimeters. Observations are in black, the BC66 model approach in red, VG80 in blue, BCVG in green, and VGc in orange. The grey shaded areas correspond to periods when the meteorological forcing is reconstructed and the blue shaded areas to the periods when data are of low quality.


With soil parameters determined in situ, the evolution of the total soil water mass is well reproduced by ISBA, whatever the model approach (Fig. 4). Skill scores (Fig. 5) are better for OPE experimental station lysimeters with a small CRMSE (below 56 kg, except for O1 with VG80). A comparison of the four model approaches shows that BC66, BCVG, and VGc exhibit better scores, with a mean CRMSE of 32, 32, and 34 kg, respectively, compared to 43 kg for VG80. VG80 also exhibits a lower R2 (0.68), compared to other experiments (0.75, 0.75, and 0.72 for BC66, BCVG, and VGc, respectively).

Figure 5Statistical scores on daily chronicles reached by BC66 (red), VG80 (blue), BCVG (green), and VGc (orange). Panel (a) shows R2 vs. CRMSE scores for the total water mass variations. Panels (c) and (d) show overall bias vs. CRMSE scores for the volumetric water content at 0.5 m depth (ω50) and 1 m depth (ω100), respectively. Panel (b) shows the overall bias vs. NASH scores for the drainage flux. All of these skill scores are presented in Appendix A. Each lysimeter is represented by its identifier. For each model approach, the ellipses represent the multivariate Student's t distribution, following Fisher (1992).


4.1.2 Water content

Time evolutions of the water content at 0.5 m are shown in Fig. 6, while the scores are presented for the two depths for which there are the most usable observations, i.e. 0.5 and 1 m (Fig. 5). For lysimeters with vegetation (G3, G4, O1, O2, and O3), the roots draw water in the summer period, which reduces the volumetric water content, causing a more pronounced contrast in the volumetric water content between winter and summer than for bare-soil lysimeters. This behaviour is well represented by ISBA because the biases are negligible (<6×10-2 m3 m−3) and dynamics are correct (R2>0.5; not shown). Still, minor differences between BC66, BCVG, and VGc appear. VG80 obtains weaker statistical scores in 65 % of the cases because soil water saturation is reached too rapidly. VGc is clearly a solution for simulating the water content when the parameter n is close to 1, by correcting the default of VG80 for these soils.

Figure 6Same as Fig. 4 but for hourly volumetric water content (m3 m−3) at 0.5 m depth.


It should be noticed that the agreement between the observed and simulated soil volumetric water content is weaker at the shallow depth (available only at OPE; 0.2 m) compared to other depths with R2<0.6 (not shown). This can be explained by the different processes which can modify the structure of the soil at the surface; the intensity of precipitation can increase the soil surface sealing (Liu et al.2011; Assouline2004), and the soil heterogeneity can increase in response to plant or biological activities (Brown et al.2000; Beven and Germann1982). Such processes are not represented in ISBA.

4.1.3 Drainage water

Drainage at 2 m depth is measured at an hourly frequency. However, to compensate for the measurement limits associated with a 0.1 mm h−1 threshold, the data are aggregated daily. The annual volume drained varies significantly between the lysimeters at the same sites (Table S2), although it is assumed that they are exposed to the same atmospheric conditions. At the GISFI experimental station, the mean annual drained water is at maximum levels for the bare-soil lysimeters G1 and G2 (>300 mm yr−1). The mean annual drained water is 2 to 3 times lower for the lysimeters covered by vegetation. At the OPE experimental station, all the lysimeters have vegetation cover, and the mean annual drained water shows a variation of only 16 %, from 300 to 363 mm yr−1. Such values are comparable to the volume drained on the bare soil of the GISFI experimental station, which is mainly due to the higher mean annual precipitation at OPE.

The daily mean annual cycles of drainage water are shown in Fig. 7. At the GISFI experimental station, drainage water occurs almost all year long for the bare-soil lysimeters G1 and G2. The well developed vegetation cover in G3 causes a decrease in both the drainage water intensity and drainage water duration. At the OPE experimental station, the drainage water occurs mainly from October to June (if the year 2016 is excluded, as a large rainfall event occurred in May–June 2016), with similar cycles for the three lysimeters. The annual cycle is well simulated, with more discrepancies for the VG80 model approach, which tends to overestimate the drainage water during some recession periods (as, for example, during spring for G1, G4, and O3).

Figure 7Same as Fig. 4 but for the daily mean annual cycles of the drainage water time series (mm d−1).


The daily drainage water is shown in Fig. 8, and the scores are given in Fig. 5. Even if the annual volume drained is higher for the OPE lysimeters, the maximum drainage water intensities over the observed period are similar at both sites; they vary between 27.4 and 34.0 mm d−1 at the OPE experimental station and between 24.0 and 33.0 mm d−1 at the GISFI site. The four model approaches reproduce the daily drainage water with relatively low biases (<0.7 mm d−1), with the worst biases being obtained by VG80, especially for the GISFI lysimeters, and confirming the results shown in Fig. 7. Dynamics are also well reproduced, as shown by the NASH scores. These scores are similar for BC66, BCVG, and VGc, with an average NASH of 0.42, 0.42, and 0.43, respectively, but with a slightly lower score, 0.35, for VG80. The worse performance of VG80 in reproducing the daily drainage water is especially true for the GISFI lysimeters.

Figure 8Same as Fig. 4 but for daily drainage water (mm d−1).


4.2 Water budget

Lysimeters give access to an estimate of the actual evapotranspiration (Schrader et al.2013; Gebler et al.2015). By neglecting the lateral runoff (not present in these lysimeters), the following simple water balance equation allows us to estimate annual evapotranspiration (E) from annual precipitation (P), annual drainage (Qdrain), and a water mass variation that is negligible (ΔW) for each lysimeter:

(7) E = P - Q drain - Δ W .

Figure 9 presents the water budget observed and simulated over the entire period for all lysimeters. The ratios of the total evapotranspiration and the drainage water to the total precipitation are expressed in percentages. At the GISFI experimental station, between 50 % and 80 % of the annual rainfall is evapotranspired, with maximal uptake for lysimeter G3, which has the densest vegetation cover. At the OPE experimental station, evapotranspiration corresponds to nearly 50 % of the rainfall.

Figure 9Water budget partition of the precipitation into drainage water and evapotranspiration (expressed in %) for lysimeters at GISFI (G1, G2, G3, and G4) and OPE (O1, O2, and O3) for each model approach. The ratio of observed drainage water to precipitation is represented by the dashed black line. Simulated drainage water is represented by colour bars for each model approach, while the simulated evapotranspiration (E) is shown in grey. The BC66 simulated drainage water is in red, VG80 in blue, BCVG in green, and VGc in orange. Drainage water simulated by additional model approaches with a homogeneous soil profile (BC66HOM and VGcHOM), with parameters estimated from usual PTFs (BC66PTF and VGcPTF), and with parameters estimated from the usual PTFs, except for b (B66PTFb) or n (VGcPTFn), estimated in situ are also shown. These last model approaches are defined in Sect. 5.


The annual budgets simulated by ISBA are rather close to the observations, but the drainage water is generally overestimated and thus the evapotranspiration underestimated on all lysimeters. The absolute difference averaged over all the lysimeters is not significant different between experiments. There are lower values for BC66 (8.1 percentage points) compared to the other model approaches (10.8 percentage points for BCVG, 11.7 for VGc, and 12.9 for VG80). This seems to underline that using the closed‐form equations from van Genuchten (1980) could favour drainage water at the expense of evapotranspiration, compared to Brooks and Corey (1966), at least over our lysimeters, and with the ISBA LSM.

4.3 Intense drainage events

In the previous sections (Sects. 4.1.3 and 4.2), drainage water has been analysed on complete chronicles, where strong daily drainage water events were detected. In order to check the ability of the four model approaches to reproduce a strong soil water dynamic, our focus is now on intense drainage water events. All daily drainage water fluxes larger than or equal to the 99th percentile of the daily drainage water distribution over the entire period are selected for each lysimeter. These Q99 values are higher for OPE experimental station lysimeters (>13 mm d−1) than for GISFI experimental station lysimeters (5.4 to 9 mm d−1). A total of 110 events on the set of the seven lysimeters is selected. In total, 75 % of these intense drainage water events appear from October to March (i.e. during the wet period when the soil is near saturation). The remainder occurs in May and June, which are associated with intense precipitation events, as generally observed in this region.

Figure 10Daily precipitation (mm d−1), hourly effective wetting saturation profile (%) observed (OBS) and simulated by BC66, VG80, BCVG, and VGc, and daily drainage water (mm d−1) observed (in black), and simulated by BC66 in red, VG80 in blue, BCVG in green, and VGc in orange during intense drainage water event in February 2016 for lysimeters O1 and G2. The NASH scores for each simulated drainage water event are also given.


Figure 10 presents the winter intense drainage water events in February 2016 for two contrasting lysimeters (O1 with vegetation and G2 with bare soil). February 2016 corresponds to a period of approximately 1 month (31 d for O1 and 27 d for G2) with some daily precipitation above 20 mm d−1 and an initially wet soil. Figure 11 presents the late-spring events in June 2016 that led to a flood event of the Seine river (Philip et al.2018). This event is characterized by large accumulated precipitation (166 mm at OPE and 210 mm at GISFI in 10 d) and intense daily precipitation, with a maximum on 30 May that reached above 40 and 70 mm d−1 at the OPE and GISFI experimental stations, respectively. Figures 10 and 11 show observed daily precipitation, hourly observed and simulated soil profile saturation, and daily observed and simulated drainage water. VG80 tends to simulate conditions that are too wet over the entire soil profile, compared to the observations. For these lysimeters, BC66, BCVG and VGc reproduce soil moisture profiles well. The drainage water simulated with VG80 is further from the observations than the other model approaches. In one case (G2 lysimeter), this VG80-simulated drainage water is occurring too early, whatever the season, while in the other case (O1 lysimeter), its dynamic is too smooth during winter. Conversely, the dynamics of these events are well reproduced in phase and maximum intensity by BC66, BCVG, and VGc, as highlighted by the good NASH scores. When the same comparison is made with the 110 selected intense drainage water events, then the scores are significantly better for BCVG and VGc. They exhibit the lowest bias (1.11 and 0.8 mm d−1) compared to the other model approaches (1.3 mm d−1 for BC66 and 1.26 mm d−1 for VG80), in addition to the highest NASH criteria (0.70 and 0.80 for BCVG and VGc, compared to 0.70 and 0.56 for BC66 and VG80).

Figure 11Same as Fig. 10 but for June 2016.


4.4 Synthesis

To summarize the results, Taylor diagrams (Taylor2001) are used to quantify the degree of correspondence between the modelled and observed behaviour in terms of three statistics, namely the Pearson correlation coefficient, the centred root mean square error, and the normalized standard deviation (Fig. 12). These scores are computed using all seven lysimeter time series, with a single result for BC66, VG80, BCVG, and VGc.

Figure 12Taylor diagrams for hourly total water mass (a), hourly volumetric water content at 0.5 m depth (b), daily drainage water (c), and daily intense drainage water events (d). Model approach BC66 is plotted in red, VG80 in blue, BCVG in green, and VGc in orange. Additional model approaches with a homogeneous soil profile (BCVGHOM and VGcHOM) are represented as open circles, with parameters estimated from the usual PTFs (BCVGPTF and VGcPTF) by pluses (+) and with parameters estimated with PTFs, except or b (BCVGPTFb) or n (VGcPTFn), which are estimated in situ, by crosses (×). The Pearson correlation coefficient, the centred root mean square error (CRMSE), and the normalized standard deviation are summarized in this diagram (Taylor2001).


For the water mass and volumetric water content (Fig. 12a and b), scores are calculated at an hourly time step, while a daily time step is used for the drainage water over both the full period and the 110 intense drainage water events (Fig. 12c and d). BC66, BCVG, and VGc obtain good results, especially for predicting water mass, volumetric water content, and drainage water during intense events. Consistent with the previous results, VG80 obtains significantly lower scores, whatever the observable. VGc exhibits a larger score in terms of the intense drainage events and the same scores as BC66 and BCVG for other variables. These results highlight that the VGc model of Iden et al. (2015) is a very interesting alternative to the VG80 model for hydrological applications with an LSM, while maintaining an approach that is integrally based on closed‐form equations from van Genuchten (1980).

5 Sensitivity model approaches

5.1 Homogeneous soil profile

As LSMs are used on regional to global scales, their soil parameters are usually derived from soil maps that generally consider a homogeneous soil profile. To evaluate the influence of the variation in the depth of soil hydrodynamic parameters on our simulations, we performed a sensitivity test with a uniform soil profile using the BC66 model (BC66HOM) and the VGc model (VGcHOM). This uniform profile is fed with the vertical mean value of each parameter for each lysimeter (crosses in Fig. 3). Using a homogeneous profile in BC66HOM and VGcHOM significantly degrades the scores in terms of the water content (Fig. 12) but has a limited impact on the simulated water budget compared to BC66 and VGc (Fig. 9). It has a stronger impact on intense drainage events than on the whole time series (Fig. 12 and Table 2). Figure 13 shows very clearly that BC66HOM fails to capture the observed soil moisture profile (drier surface conditions and wetter with depth), which is in contrast with BC66 (Figs. 10 and 11). VGcHOM exhibits the same behaviour (Fig. S5). The drainage water dynamics are less impacted during winter than during spring, as shown by the BC66HOM NASH scores, compared to BC66 (and accordingly for VGcHOM and VGc). Indeed, lysimeter soils simulated with a homogeneous profile appear wetter than the observations and reference simulations, especially in spring. This wetter state thus induces the reactivity of the drainage water during this period that is too intense, thereby reducing the simulated skill scores compared to the reference simulations. This fact underlines that using a homogeneous soil profile sometimes fails to correctly simulate the drainage dynamic in the studied lysimeters for some conditions.

Figure 13Same as Figs. 10 and 11 but for sensitivity model approaches with a homogeneous soil profile (BC66HOM), with soil parameters from the usual PTFs (BC66PTF), and with parameters from the usual PTFs, except for b, estimated in situ (BC66PTFb).


Table 2NASH scores for the simulated drainage water (Qdrain) over the seven lysimeters and the entire period and during intense drainage water events (Qint). Model approaches are shown with the soil hydrodynamic parameters set, with a homogeneous vertical profile (BC66HOM and VGcHOM), computed using the usual PTFs (BC66PTF and VGcPTF), or derived from observations (BC66 and VGc). The NASH of the additional model approaches with one parameter (n, ψsat, ωsat, or ksat) that keeps the reference values is also given.

NA: not available

Download Print Version | Download XLSX

In order to determine which hydrodynamical parameter (via its vertical profile) has the largest impact on the simulations, we performed additional sensitivity tests with homogeneous soil profiles for all parameters, except for one that keeps its estimated heterogeneous profile. These tests are performed with all parameters being homogeneous, except one, with either ωsat, b, or ψsat for BC66 and ωsat, n, or α for VGc.

Using the seven lysimeters, we are able to complete the drainage water time series and all the selected intense drainage water events, as in Sect. 4.4. Table 2 shows that b and n are the most important parameters, as accounting for their heterogeneous profiles improves the NASH score more, compared to the other parameters. Their NASH scores are, in addition, very close to the BC66 and VGc references. These tests demonstrate the importance of b (and therefore n) for accurately simulating the drainage water dynamic and intense drainage water events. This finding is in agreement with previous studies (Ritter et al.2003) that demonstrated a strong sensitivity of the simulated drainage water to n (and therefore b) and a lower sensitivity to ksat.

5.2 Usual pedotransfer functions

LSMs commonly use PTFs to derive soil hydrodynamic parameters from soil textural information. As shown in Fig. 3, the soil parameters estimated from our measurements can be very different from those derived from six usual PTFs (Clapp and Hornberger1978; Cosby et al.1984; Carsel and Parrish1988; Wösten and van Genuchten1988; Vereecken et al.1989; Weynants et al.2009). This is especially true for the b and n parameters. To investigate the impacts of such differences, sensitivity tests noted that BC66PTF and VGcPTF are performed so that the soil hydrodynamic parameters in each soil horizon are derived from the mean of these PTFs. Impacts on the simulated water budget are clear but not obvious enough to comment on (Fig. 9), even if BC66PTF and VGcPTF tend to increase the drainage water at the expense of evapotranspiration compared to BC66 and VGc in the GISFI lysimeters, while the opposite behaviour (lower drainage water and larger evapotranspiration) is found for the OPE lysimeters. All skill scores are drastically degraded in terms of water mass variations, volumetric water content at 0.5 m depth, daily drainage water, and intense drainage events (Fig. 12 and Table 2). These weaknesses are highlighted over the selected February and June 2016 events (Figs. 13 and S5). The soil water profile simulated by BC66PTF and VGcPTF is strongly underestimated compared to the observations, especially for the G2 lysimeter. This weakness induces a significant delay in the simulated drainage water, compared to observations and to other simulations.

Once again, the b and n parameters seem to be the key to this weakness. Indeed, we performed a set of tests using the BC66PTF and VGcPTF configuration, except for one parameter, either b (or n), ωsat, ψsat (or α), or ksat, that keeps its estimated in situ value. The use of the in situ estimated values of b and n tends to balance the water budget toward its partitioning of the reference (Fig. 9), especially with VGc. It drastically improves all skill scores compared to BC66PTF and VGcPTF (Fig. 12). This is also the case for the simulation of the soil water profile and especially the soil drainage water dynamic (Table 2; Figs. 13 and S5).

5.3 Lower boundary conditions

LSMs commonly use free-drainage lower boundary conditions (LBCs) on field sites or in regional- and global-scale application. We wondered if such an LBC is able to impact our results, compared to the seepage face LBC classically used to simulate lysimeters. Four simulations were performed using the same four model approaches as previously but with a free-drainage LBC instead of a seepage face LBC.

Figure 14Mean annual cycles of the hourly effective wetting saturation profile (%) and daily drainage water (mm d−1) observed (OBS in black) and simulated by BC66 (left panels) and VGc (right panels), with a seepage face LBC (BC66SEEPAGE in red and VGcSEEPAGE in orange) and a free-drainage LBC (BC66FREE and VGcFREE in grey) over lysimeters O1 and G2. The NASH scores for each simulated drainage water are also given.


Figure 14 shows the mean daily annual cycles of the soil moisture profiles and the drainage water observed and simulated by BC66 and VGc, with the two LBC approaches (BC66FREE and VGcFREE compared to BC66SEEPAGE and VGcSEEPAGE) for the O1 and G2 lysimeters. Soil moisture profiles are sensibly the same, even if, logically, the deep layers are slightly drier because the deepest layer is always saturated with a seepage face LBC. This fact is confirmed for all model approaches (BC66, VG80, BCVG, and VGc) which exhibit the same skill score in terms of water mass variations and volumetric water content at 0.5 m, whatever the LBC (comparing Figs. 12 and S6). The simulated total water budgets are also relatively unchanged (Fig. S7).

In terms of the daily drainage, the simulated response is not very different, whatever the LBC (Fig. 14). To use a seepage face LBC seems effectively more adequate for simulating drainage water in our lysimeters during recession periods, especially with BC66 (that explains the larger NASH for BC66SEEPAGE compared to BC66FREE for the G2 lysimeter). Some peaks of drainage water are also more pronounced with such an LBC. Comparing the Taylor diagrams (Figs. 12 and S6) underlines that all model approaches (BC66, VG80, BCVG, and VGc) reproduce daily drainage water and intense drainage events with the same accuracy, whatever the LBC.

6 Discussion and conclusion

This study uses time series (up to 7 years) from several lysimeters to evaluate the dynamics of water transfer in the unsaturated zone simulated with the land surface model ISBA. These observations allow the derivation of the heterogeneous profile of soil parameters. Although the original version of ISBA performed well, a set of four water closure relationships, which estimate the evolution of soil hydrodynamic properties with soil moisture, are tested. A comparison of these four relationships shows that, when soil parameters and meteorologic forcing are known, ISBA reproduces the evolution of soil hydrology and vegetation processes reasonably well. The choice to use a seepage face LBC in ISBA (as commonly done to simulate lysimeters) instead of its usual free-drainage LBC has a very little impact on the results.

The simulation using the VG80 water closure relationships exhibits more difficulty when reproducing the soil water profile and the drainage water dynamic, in particular during intense drainage water events, than the original ISBA version using the BC66 equations. This is partly linked to the limitation of the VG80 hydraulic conductivity function for n close to 1. The BCVG model approach combines the soil matrix potential function from VG80 with the hydraulic conductivity function from BC66, and the VGc model approach of Iden et al. (2015) that includes a maximum pore radius on the VG80 hydraulic conductivity curve solves these problems and even seems to be able to improve the simulation of the soil drainage water dynamic compared to BC66. However, the l parameter in the VG80 hydraulic conductivity (Eq. 3), and thus VGc (Eq. 5), is difficult to set, even with direct observations. In our study, it is fixed at the classical 0.5 value for some lysimeters but can vary drastically in others, which is consistent with the literature (Wösten and van Genuchten1988; Wösten et al.2001; Schaap et al.2001), underlying the difficulty with respect to estimating this parameter for regional- to global-scale applications.

The observations show that the soil hydrodynamic parameters in each lysimeter are strongly heterogeneous with depth, while LSMs generally use homogeneous profiles. Using additional sensitivity tests with such homogeneous profiles, we found that, even if the simulated soil water and drainage water dynamics remain acceptable compared to the observations, then all of the skill scores are worsened (especially for the soil water profile) compared to the model approaches with a heterogeneous profile. This finding supports the need to account for the vertical heterogeneity of soil hydrodynamic parameters (King et al.1999; Mirus2015; Hengl et al.2017; Vogel2019; Fatichi et al.2020; Bauser et al.2020; Gebler et al.2017) to improve the simulation of soil water and drainage water dynamics (Stieglitz et al.1997; Mohanty and Zhu2007; Decharme et al.2011; Vereecken et al.2019). This is a challenge for simulating groundwater recharge on regional and global scales.

We also found that parameters b and n, which characterize the shape of the soil water retention function, derived from the observations differ significantly from those derived from PTFs commonly used in LSMs (Clapp and Hornberger1978; Cosby et al.1984; Carsel and Parrish1988; Wösten and van Genuchten1988; Vereecken et al.1989; Weynants et al.2009). Sensitivity tests show that the values of b and n derived from usual PTFs are not suitable for simulating the drainage water dynamic, at least for the seven lysimeters used in this study. In addition, these parameters exhibit the largest heterogeneity with soil depth. Neglecting this behaviour contributes to degrading the simulated drainage water dynamic. Note that this heterogeneous behaviour of b and n is still under consideration in the literature. As in our study, some authors observed a decrease in n with soil depth (Ritter et al.2003; Jhorar et al.2004; Schwärzel et al.2006), while some others showed an increase (Groh et al.2018) or even no change (Schneider et al.2021). In any case, b and n could be the key parameters for correctly simulating the drainage water dynamic and groundwater recharge with LSMs. We recognize that this last assumption has to be confirmed in many other experimental field sites. Indeed, this study is based on two experimental sites with similar climatic conditions and with low-intensity precipitation events compared to other regions. It would be interesting to conduct additional studies in other contrasting climates.

In the context of LSMs that can be used at the regional or global scale, the major challenge is to simplify the numerical simulations and parameter calibration to be applicable in different contexts, without specific calibration, and to reproduce the water cycle as much as possible. This study at the local scale increases the confidence that LSMs are powerful tools for simulating the recharge of groundwater in different environmental conditions, with many soils and vegetation cover types, and can therefore be used for many applications in hydrology at both the regional and global scales. The sensitivity of our LSM to the soil heterogeneity or the value of some hydrodynamical parameters underlines, however, that this remains a challenge, even with the advent of global databases describing the vertical profile of the soil properties at depths greater than 1 m (Poggio et al.2021).

Appendix A: Statistical scores

In this study, we used the following skill scores considering the simulated, Si, and the observed, Oi, data defined at N discrete points in time, including:

  • the overall bias,

    (A1) Bias = S - O with S = 1 N i = 1 N S i and O = 1 N i = 1 N O i ;
  • the centred root mean square error,

    (A2) CRMSE = 1 N i = 1 N S i - S - O i - O 2 ;
  • the Pearson correlation coefficient (or the coefficient of determination),

    (A3) R = 1 N i = 1 N S i - S - O i - O σ s σ o ,

    where σs and σo are the standard deviations of Si and Oi, respectively. The coefficient of determination, R2, is the square of the Pearson correlation; and

  • the Nash and Sutcliffe (1970) efficiency criterion,

    (A4) NASH = 1 - i = 1 N S i - O i 2 i = 1 N O i - O 2 .
Code availability

The SURFEX v8.1 code, including the ISBA code, used in this study is freely available at (Le Moigne, 2018).

Data availability

The lysimetric data are not publicly available and are provided only upon request to the corresponding author (to be checked with the data owners Noële Enjelvin, for GISFI, and Paul-Olivier Redon, for OPE).


The supplement related to this article is available online at:

Author contributions

The article was written by AS, with contributions from all co-authors. AS, BD, and FH developed the study design. NE, POR, and PFC collected the lysimeter data and assisted in their use. AS and BD developed new parameterizations. AS analysed the data and ran the models, with general support from BD and FH. Specific support was provided by CD for parameterizing vegetation and by PLM for computing the weather forcing.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors are very grateful to Corinne Leyval from LIEC Nancy and Geoffroy Séré from LSE for all the useful discussions and their willingness to provide access to the datasets from GISFI. The authors thank GISFI for providing field site data at the experimental station in Homécourt and in particular the technical staff, Cyrielle Boone. The OPE dataset was provided by Andra, with the help of Catherine Galy from Andra. The authors are also very grateful to Wolfgang Durner, who suggested the test of the VGc approach. The authors extend their thanks to IR OZCAR, whose infrastructure enabled the connection between research from different topics at the various sites.

Financial support

Antoine Sobaga has been funded by the Centre National de Recherches Météorologiques of Météo-France, the Agence de l’Eau Seine Normandie, and the Laboratoire de Géologie de l’Ecole Normale Supérieure of Paris.

Review statement

This paper was edited by Philippe Ackerer and reviewed by two anonymous referees.


Alcamo, J., Döll, P., Henrichs, T., Kaspar, F., Lehner, B., Rösch, T., and Siebert, S.: Development and testing of the WaterGAP 2 global model of water use and availability, Hydrolog. Sci. J., 48, 317–337, 2003. a

Assouline, S.: Rainfall-Induced Soil Surface Sealing: A Critical Review of Observations, Conceptual Models, and Solutions, Vadose Zone J., 3, 570–591,, 2004. a

Bauser, H., Riedel, L., Berg, D., Troch, P., and Roth, K.: Challenges with effective representations of heterogeneity in soil hydrology based on local water content measurements, Vadose Zone J., 19, e20040,, 2020. a

Beven, K. and Germann, P.: Macropores and Water Flow in Soils, Water Resour. Res., 18, 1311–1325,, 1982. a

Blöschl, G., Bierkens, M. F., Chambel, A., et al.: Twenty-three unsolved problems in hydrology (UPH) – a community perspective, Hydrolog. Sci. J., 64, 1141–1158, 2019. a

Blyth, E., Gash, J., Lloyd, A., Pryor, M., Weedon, G. P., and Shuttleworth, J.: Evaluating the JULES land surface model energy fluxes using FLUXNET data, J. Hydrometeorol., 11, 509–519, 2010. a

Boone, A. and Etchevers, P.: An Intercomparison of Three Snow Schemes of Varying Complexity Coupled to the Same Land Surface Model: Local-Scale Evaluation at an Alpine Site, J. Hydrometeorol., 2, 374–394,<0374:AIOTSS>2.0.CO;2, 2001. a

Boone, A., Masson, V., Meyers, T., and Noilhan, J.: The Influence of the Inclusion of Soil Freezing on Simulations by a Soil–Vegetation–Atmosphere Transfer Scheme, J. Appl. Meteorol., 39, 1544–1569,<1544:TIOTIO>2.0.CO;2, 2000. a, b, c, d, e, f

Braud, I., Dantas-Antonino, A., Vauclin, M., Thony, J., and Ruelle, P.: A simple soil-plant-atmosphere transfer model (SiSPAT) development and field verification, J. Hydrol., 166, 213–250,, 1995. a, b, c, d

Braud, I., Chaffard, V., Coussot, C., Galle, S., Juen, P., Alexandre, H., Baillion, P., Battais, A., Boudevillain, B., Branger, F., Brissebrat, G., Cailletaud, R., Cochonneau, G., Decoupes, R., Desconnets, J.-C., Dubreuil, A., Fabre, J., Gabillard, S., Gerard, M.-F., Grellet, S., Herrmann, A., Laarman, O., Lajeunesse, E., Le Henaff, G., Lobry, O., Mauclerc, A., Paroissien, J.-B., Pierret, M.-C., Silvera, N., and Squividant, H.: Building the information system of the French Critical Zone Observatories network: Theia/OZCAR-IS, Hydrolog. Sci. J., 67, 1–19,, 2020. a

Bredehoeft, J. D.: The water budget myth revisited: why hydrogeologists model, Groundwater, 40, 340–345,, 2002. a

Brooks, R. H. and Corey, A. T.: Properties of Porous Media Affecting Fluid Flow, J. Irrig. Drain. Div., 92, 61–88,, 1966. a, b, c, d, e, f, g

Brown, G. G., Barois, I., and Lavelle, P.: Regulation of soil organic matter dynamics and microbial activityin the drilosphere and the role of interactionswith other edaphic functional domains Paper presented at the 16th World Congress of Soil Science, 20–26 August 1998, Montpellier, France, Eur. J. Soil Biol., 36, 177–198,, 2000. a

Brunner, P., Therrien, R., Renard, P., Simmons, C. T., and Franssen, H.-J. H.: Advances in understanding river-groundwater interactions, Rev. Geophys., 55, 818–854,, 2017. a

Burdine, N.: Relative Permeability Calculations From Pore Size Distribution Data, J. Petrol. Technol., 5, 71–78,, 1953. a, b

Calvet, J.-C.: Investigating soil and atmospheric plant water stress using physiological and micrometeorological data, Agr. Forest Meteorol., 103, 229–247,, 2000. a

Calvet, J.-C. and Soussana, J.-F.: Modelling CO2-enrichment effects using an interactive vegetation SVAT scheme, Agr. Forest Meteorol., 108, 129–152,, 2001. a

Calvet, J.-C., Noilhan, J., Roujean, J.-L., Bessemoulin, P., Cabelguenne, M., Olioso, A., and Wigneron, J.-P.: An interactive vegetation SVAT model tested against data from six contrasting sites, Agr. Forest Meteorol., 92, 73–95,, 1998. a, b

Carrick, S., Rogers, G., Cameron, K., Malcolm, B., and Payne, J.: Testing large area lysimeter designs to measure leaching under multiple urine patches, NZ J. Agricult. Res., 60, 205–215, 2017. a

Carsel, R. F. and Parrish, R. S.: Developing joint probability distributions of soil water retention characteristics, Water Res. Res., 24, 755–769,, 1988. a, b, c, d, e

Clapp, R. B. and Hornberger, G. M.: Empirical equations for some soil hydraulic properties, Water Resour. Res., 14, 601–604,, 1978. a, b, c, d

Cosby, B. J., Hornberger, G. M., Clapp, R. B., and Ginn, T. R.: A Statistical Exploration of the Relationships of Soil Moisture Characteristics to the Physical Properties of Soils, Water Resour. Res., 20, 682–690,, 1984. a, b, c, d

Darvishzadeh, R., Skidmore, A., Schlerf, M., Atzberger, C., Corsi, F., and Cho, M.: LAI and chlorophyll estimation for a heterogeneous grassland using hyperspectral measurements, ISPRS J. Photogram. Remote Sens., 63, 409–426,, 2008. a

Deardorff, J. W.: Efficient prediction of ground surface temperature and moisture, with inclusion of a layer of vegetation, J. Geophys. Res.-Oceans, 83, 1889–1903, 1978. a

Decharme, B., Boone, A., Delire, C., and Noilhan, J.: Local evaluation of the Interaction between Soil Biosphere Atmosphere soil multilayer diffusion scheme using four pedotransfer functions, J. Geophys. Res., 116, D20126,, 2011. a, b, c, d, e, f, g

Decharme, B., Brun, E., Boone, A., Delire, C., Le Moigne, P., and Morin, S.: Impacts of snow and organic soils parameterization on northern Eurasian soil temperature profiles simulated by the ISBA land surface model, The Cryosphere, 10, 853–877,, 2016. a, b

Decharme, B., Delire, C., Minvielle, M., Colin, J., Vergnes, J.-P., Alias, A., Saint-Martin, D., Séférian, R., Sénési, S., and Voldoire, A.: Recent Changes in the ISBA-CTRIP Land Surface System for Use in the CNRM-CM6 Climate Model and in Global Off-Line Hydrological Applications, J. Adv. Model. Earth Syst., 11, 1207–1252,, 2019. a, b, c

Delire, C., Séférian, R., Decharme, B., Alkama, R., Calvet, J.-C., Carrer, D., Gibelin, A.-L., Joetzjer, E., Morel, X., Rocher, M., and Tzanos, D.: The Global Land Carbon Cycle Simulated With ISBA-CTRIP: Improvements Over the Last Decade, J. Adv. Model. Earth Syst., 12, e2019MS001886,, 2020. a

De Rosnay, P., Polcher, J., Laval, K., and Sabre, M.: Integrated parameterization of irrigation in the land surface model ORCHIDEE. Validation over Indian Peninsula, Geophys. Res. Lett., 30, 1986,, 2003. a

Dettmann, U., Andrews, F., and Donckels, B.: Package “SoilHyP”, (last access: 1 June 2023), 2022. a

Döll, P. and Fiedler, K.: Global-scale modeling of groundwater recharge, Hydrol. Earth Syst. Sci., 12, 863–885,, 2008. a

Dourado Neto, D., de Jong van Lier, Q., van Genuchten, M., Reichardt, K., Metselaar, K., and Nielsen, D.: Alternative Analytical Expressions for the General van Genuchten–Mualem and van Genuchten–Burdine Hydraulic Conductivity Models, Vadose Zone J., 10, 618–623,, 2011. a

Durner, W.: Hydraulic conductivity estimation for soils with heterogeneous pore structure, Water Resour. Res., 30, 211–223, 1994. a

Farthing, M. W. and Ogden, F. L.: Numerical Solution of Richards' Equation: A Review of Advances and Challenges, Soil Sci. Soc. Am. J., 81, 1257–1269,, 2017. a

Fatichi, S., Or, D., Walko, R., Vereecken, H., Young, M., Ghezzehei, T., Hengl, T., Kollet, S., Agam, N., and Avissar, R.: Soil structure is an important omission in Earth System Models, Nat. Commun., 11, 522,, 2020. a

Fisher, R. A.: Statistical methods for research workers, in: Breakthroughs in statistics, Springer, 66–70, ISBN: 9788130701332, 1992. a

Fuentes, C., Haverkamp, R., and Parlange, J.-Y.: Parameter constraints on closed-form soilwater relationships, J. Hydrol., 134, 117–142,, 1992. a, b, c, d, e

Gaillardet, J., Braud, I., Hankard, F., Anquetin, S., Bour, O., Dorfliger, N., de Dreuzy, J., Galle, S., Galy, C., Gogo, S., Gourcy, L., Habets, F., Laggoun, F., Longuevergne, L., Le Borgne, T., Naaim-Bouvet, F., Nord, G., Simonneaux, V., Six, D., Tallec, T., Valentin, C., Abril, G., Allemand, P., Arènes, A., Arfib, B., Arnaud, L., Arnaud, N., Arnaud, P., Audry, S., Comte, V. B., Batiot, C., Battais, A., Bellot, H., Bernard, E., Bertrand, C., Bessière, H., Binet, S., Bodin, J., Bodin, X., Boithias, L., Bouchez, J., Boudevillain, B., Moussa, I. B., Branger, F., Braun, J. J., Brunet, P., Caceres, B., Calmels, D., Cappelaere, B., Celle-Jeanton, H., Chabaux, F., Chalikakis, K., Champollion, C., Copard, Y., Cotel, C., Davy, P., Deline, P., Delrieu, G., Demarty, J., Dessert, C., Dumont, M., Emblanch, C., Ezzahar, J., Estèves, M., Favier, V., Faucheux, M., Filizola, N., Flammarion, P., Floury, P., Fovet, O., Fournier, M., Francez, A. J., Gandois, L., Gascuel, C., Gayer, E., Genthon, C., Gérard, M. F., Gilbert, D., Gouttevin, I., Grippa, M., Gruau, G., Jardani, A., Jeanneau, L., Join, J. L., Jourde, H., Karbou, F., Labat, D., Lagadeuc, Y., Lajeunesse, E., Lastennet, R., Lavado, W., Lawin, E., Lebel, T., Le Bouteiller, C., Legout, C., Lejeune, Y., Le Meur, E., Le Moigne, N., Lions, J., Lucas, A., Malet, J. P., Marais-Sicre, C., Maréchal, J. C., Marlin, C., Martin, P., Martins, J., Martinez, J. M., Massei, N., Mauclerc, A., Mazzilli, N., Molénat, J., Moreira-Turcq, P., Mougin, E., Morin, S., Ngoupayou, J. N., Panthou, G., Peugeot, C., Picard, G., Pierret, M. C., Porel, G., Probst, A., Probst, J. L., Rabatel, A., Raclot, D., Ravanel, L., Rejiba, F., René, P., Ribolzi, O., Riotte, J., Rivière, A., Robain, H., Ruiz, L., Sanchez-Perez, J. M., Santini, W., Sauvage, S., Schoeneich, P., Seidel, J. L., Sekhar, M., Sengtaheuanghoung, O., Silvera, N., Steinmann, M., Soruco, A., Tallec, G., Thibert, E., Lao, D. V., Vincent, C., Viville, D., Wagnon, P., and Zitouna, R.: OZCAR: The French Network of Critical Zone Observatories, Vadose Zone J., 17, 180067,, 2018. a

Garrigues, S., Olioso, A., Carrer, D., Decharme, B., Calvet, J., Martin, E., Moulin, S., and Marloie, O.: Impact of climate, vegetation, soil and crop management variables on multi-year ISBA-A-gs simulations of evapotranspiration over a mediterranean crop site, Geosci. Model Dev., 8, 3033–3053,, 2015. a

Garrigues, S., Boone, A. A., Decharme, B., Olioso, A., Albergel, C., Calvet, J.-C., Moulin, S., Buis, S., and Martin, E.: Impacts of the Soil Water Transfer Parameterization on the Simulation of Evapotranspiration over a 14-Year Mediterranean Crop Succession, J. Hydrometeorol., 19, 3–25,, 2018. a

Gebler, S., Hendricks Franssen, H. J., Pütz, T., Post, H., Schmidt, M., and Vereecken, H.: Actual evapotranspiration and precipitation measured by lysimeters: A comparison with eddy covariance and tipping bucket, Hydrol. Earth Syst. Sci., 19, 2145–2161,, 2015. a

Gebler, S., Hendricks Franssen, H. J., Kollet, S. J., Qu, W., and Vereecken, H.: High resolution modelling of soil moisture patterns with TerrSysMP: A comparison with sensor network data, J. Hydrol., 547, 309–331,, 2017. a

Gelati, E., Decharme, B., Calvet, J.-C., Minvielle, M., Polcher, J., Fairbairn, D., and Weedon, G. P.: Hydrological assessment of atmospheric forcing uncertainty in the Euro-Mediterranean area using a land surface model, Hydrol. Earth Syst. Sci., 22, 2091–2115,, 2018. a

Gibelin, A.-L., Calvet, J.-C., Roujean, J.-L., Jarlan, L., and Los, S. O.: Ability of the land surface model ISBA-A-gs to simulate leaf area index at the global scale: Comparison with satellites products, J. Geophys. Res.-Atmos., 111, D18102,, 2006. a

Goudriaan, J., van Laar, H. H., van Keulen, H., and Louwerse, W.: Photosynthesis, CO2 and Plant Production, Springer US, Boston, MA, 107–122,, 1985. a

Groh, J., Stumpp, C., Lücke, A., Pütz, T., Vanderborght, J., and Vereecken, H.: Inverse Estimation of Soil Hydraulic and Transport Parameters of Layered Soils from Water Stable Isotope and Lysimeter Data, Vadose Zone J., 17, 170168,, 2018. a

Haddeland, I., Clark, D. B., Franssen, W., Ludwig, F., Voß, F., Arnell, N. W., Bertrand, N., Best, M., Folwell, S., Gerten, D., Gomes, S., Gosling, S. N., Hagemann, S., Hanasaki, N., Harding, R., Heinke, J., Kabat, P., Koirala, S., Oki, T., Polcher, J., Stacke, T., Viterbo, P., Weedon, G. P., and Yeh, P.: Multimodel Estimate of the Global Terrestrial Water Balance: Setup and First Results, J. Hydrometeorol., 12, 869–884,, 2011. a

Haines, W. B.: Studies in the physical properties of soil. V. The hysteresis effect in capillary properties, and the modes of moisture distribution associated therewith, J. Agricult. Sci., 20, 97–116, 1930. a

Harbaugh, A. W.: MODFLOW-2005, the US Geological Survey modular ground-water model: the ground-water flow process, US Department of the Interior, US Geological Survey, Reston, VA, (last access: 1 June 2023), 2005. a

Hazelton, P. and Murphy, B.: Interpreting soil test results: What do all the numbers mean?, CSIRO publishing, ISBN: 9781523109418, 2016. a

Healy, R. and Cook, P.: Using Groundwater Levels to Estimate Recharge, Hydrogeol. J., 10, 91–109,, 2002. a

Hengl, T., Mendes de Jesus, J., Heuvelink, G. B. M., Ruiperez Gonzalez, M., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R. A., Batjes, N. H., Leenaars, J. G. B., Ribeiro, E., Wheeler, I., Mantel, S., and Kempen, B.: SoilGrids250m: Global gridded soil information based on machine learning, PLOS ONE, 12, 1–40,, 2017. a

Huot, H., Morel, J.-L., and Simonnot, M.-O.: Pedogenetic Trends in Soils formed in Technogenic Materials, Soil Sci., 180, 182–192,, 2015. a

Iden, S. C., Peters, A., and Durner, W.: Improving prediction of hydraulic conductivity by constraining capillary bundle models to a maximum pore size, Adv. Water Resour., 85, 86–92, 2015. a, b, c, d, e, f, g, h, i, j, k, l

Ippisch, O., Vogel, H.-J., and Bastian, P.: Validity limits for the van Genuchten–Mualem model and implications for parameter estimation and numerical simulation, Adv. Water Resour., 29, 1780–1789, 2006. a, b, c, d, e

Jackson, R. B., Canadell, J., Ehleringer, J. R., Mooney, H. A., Sala, O. E., and Schulze, E. D.: A global analysis of root distributions for terrestrial biomes, Oecologia, 108, 389–411,, 1996. a

Jacobs, C., van den Hurk, B., and de Bruin, H.: Stomatal behaviour and photosynthetic rate of unstressed grapevines in semi-arid conditions, Agr. Forest Meteorol., 80, 111–134,, 1996. a

Jhorar, R., Van Dam, J., Bastiaanssen, W., and Feddes, R.: Calibration of effective soil hydraulic parameters of heterogeneous soil profiles, J. Hydrol., 285, 233–247, 2004. a

Joetzjer, E., Delire, C., Douville, H., Ciais, P., Decharme, B., Carrer, D., Verbeeck, H., De Weirdt, M., and Bonal, D.: Improving the ISBACC land surface model simulation of water and carbon fluxes and stocks over the Amazon forest, Geosci. Model Dev., 8, 1709–1727,, 2015. a, b

Kattge, J., Bönisch, G., Díaz, S., Lavorel, S., Prentice, I. C., Leadley, P., Tautenhahn, S., and Werner, G.: TRY plant trait database – enhanced coverage and open access, Global Change Biol., 26, 119–188,, 2020. a

Keshavarzi, M., Baker, A., Kelly, B., and Andersen, M.: River–groundwater connectivity in a karst system, Wellington, New South Wales, Australia, Hydrogeol. J., 25, 557,, 2016. a

King, D. D., Stengel, P., and Jamagne, M.: Soil mapping and soil monitoring: state of progress and use in France, (last access: 1 June 2023), 1999. a

Kosugi, K.: Three-parameter lognormal distribution model for soil water retention, Water Resour. Res., 30, 891–901, 1994. a

Leij, F. J., Haverkamp, R., Fuentes, C., Zatarain, F., and Ross, P. J.: Soil Water Retention, Soil Sci. Soc. Am. J., 69, 1891–1901,, 2005. a

Lemaire, J., Mora, V., Faure, P., Hanna, K., Buès, M., and Simonnot, M.-O.: Chemical oxidation efficiency for aged, PAH-contaminated sites: An investigation of limiting factors, J. Environ. Chem. Eng., 7, 103061,, 2019. a

Le Moigne P.: SURFEX scientific documentation, SURFEX v8.1, Issue 3, Météo-France, Toulouse, France, (last access: 1 May 2023), 2018. 

Le Moigne, P., Besson, F., Martin, E., Boé, J., Boone, A., Decharme, B., Etchevers, P., Faroux, S., Habets, F., Lafaysse, M., Leroux, D., and Rousset-Regimbeau, F.: The latest improvements with SURFEX v8.0 of the Safran–Isba–Modcou hydrometeorological model for France, Geosci. Model Dev., 13, 3925–3946,, 2020. a, b, c

Lenhard, R. J., Parker, J. C., and Mishra, S.: On the Correspondence between Brooks‐Corey and van Genuchten Models, J. Irrig. Drain. Eng., 115, 744–751,, 1989. a, b, c

Leyval, C., Beguiristain, T., Enjelvin, N., Faure, P., Lorgeoux, C., and Putz, T.: Long term fate of multi-contamination within the context of natural attenuation and climatic conditions: a collaborative study in TERENO and SFI lysimeters, in: TERENO International Conference, Berlin, Germany, 8 October 2018, (last access: 1 June 2023), 2018. a

Liu, H., Lei, T. W., Zhao, J., Yuan, C. P., Fan, Y. T., and Qu, L. Q.: Effects of rainfall intensity and antecedent soil water content on soil infiltrability under rainfall conditions using the run off-on-out method, J. Hydrol., 396, 24–32,, 2011. a

Ma, Q., Hook, J. E., and Ahuja, L. R.: Influence of three-parameter conversion methods between van Genuchten and Brooks-Corey Functions on soil hydraulic properties and water-balance predictions, Water Resour. Res., 35, 2571–2578,, 1999. a

Mahfouf, J. F. and Noilhan, J.: Comparative Study of Various Formulations of Evaporations from Bare Soil Using In Situ Data, J. Appl. Meteorol. Clim., 30, 1354–1365,<1354:CSOVFO>2.0.CO;2, 1991. a

Maxwell, R. M. and Condon, L. E.: Connections between groundwater flow and transpiration partitioning, Science, 353, 377–380,, 2016. a

Meyboom, P.: Estimating ground-water recharge from stream hydrographs, J. Geophys. Res., 66, 1203–1214,, 1961. a

Mirus, B. B.: Evaluating the importance of characterizing soil structure and horizons in parameterizing a hydrologic process model, Hydrol. Process., 29, 4611–4623,, 2015. a

Moeck, C., Grech-Cumbo, N., Podgorski, J., Bretzler, A., Gurdak, J. J., Berg, M., and Schirmer, M.: A global-scale dataset of direct natural groundwater recharge rates: A review of variables, processes and relationships, Sci. Total Environ., 717, 137042,, 2020. a, b

Mohanty, B. P. and Zhu, J.: Effective Hydraulic Parameters in Horizontally and Vertically Heterogeneous Soils for Steady-State Land–Atmosphere Interaction, J. Hydrometeorol., 8, 715–729,, 2007. a

Monserie, M.-F., Watteau, F., Villemin, G., Ouvrard, S., and Morel, J.-L.: Technosol genesis: identification of organo-mineral associations in a young Technosol derived from coking plant waste materials, J. Soil. Sediment., 9, 537–546,, 2009. a

Morel, X., Decharme, B., Delire, C., Krinner, G., Lund, M., Hansen, B. U., and Mastepanov, M.: A New Process-Based Soil Methane Scheme: Evaluation Over Arctic Field Sites With the ISBA Land Surface Model, J. Adv. Model. Earth Syst., 11, 293–326,, 2019. a

Morel-Seytoux, H. J., Meyer, P. D., Nachabe, M., Tourna, J., van Genuchten, M. T., and Lenhard, R. J.: Parameter Equivalence for the Brooks-Corey and Van Genuchten Soil Characteristics: Preserving the Effective Capillary Drive, Water Resour. Res., 32, 1251–1258,, 1996. a

Mualem, Y.: A new model for predicting the hydraulic conductivity of unsaturated porous media, Water Resour. Res., 12, 513–522,, 1976. a, b, c

Muñoz-Sabater, J., Dutra, E., Agustí-Panareda, A., Albergel, C., Arduini, G., Balsamo, G., Boussetta, S., Choulga, M., Harrigan, S., Hersbach, H., Martens, B., Miralles, D. G., Piles, M., Rodríguez-Fernández, N. J., Zsoter, E., Buontempo, C., and Thépaut, J.-N.: ERA5-Land: a state-of-the-art global reanalysis dataset for land applications, Earth Syst. Sci. Data, 13, 4349–4383,, 2021. a

Nash, J. E. and Sutcliffe, J. V.: River flow forecasting through conceptual models part I – A discussion of principles, J. Hydrol., 10, 282–290, 1970. a, b

Nimmo, J. R.: Theory for Source-Responsive and Free-Surface Film Modeling of Unsaturated Flow, Vadose Zone J., 9, 295–306,, 2010. a

Noilhan, J. and Mahfouf, J. F.: The ISBA land surface parameterisation scheme, Global Planet. Change, 13, 145–159, 1996. a

Noilhan, J. and Planton, S.: A Simple Parameterization of Land Surface Processes for Meteorological Models, Mon. Weather Rev., 117, 536–549,<0536:ASPOLS>2.0.CO;2, 1989. a, b, c

Ouvrard, S., Barnier, C., Bauda, P., Beguiristain, T., Biache, C., Bonnard, M., Caupert, C., Cébron, A., Cortet, J., Cotelle, S., Dazy, M., Faure, P., Masfaraud, J. F., Nahmani, J., Palais, F., Poupin, P., Raoult, N., Vasseur, P., Morel, J. L., and Leyval, C.: In Situ Assessment of Phytotechnologies for Multicontaminated Soil Management, Int. J. Phytoremed., 13, 245–263,, 2011. a, b

Philip, S., Kew, S. F., van Oldenborgh, G. J., Aalbers, E., Vautard, R., Otto, F., Haustein, K., Habets, F., and Singh, R.: Validation of a rapid attribution of the May/June 2016 flood-inducing precipitation in France to climate change, J. Hydrometeorol., 19, 1881–1898,, 2018. a

Poggio, L., de Sousa, L. M., Batjes, N. H., Heuvelink, G. B. M., Kempen, B., Ribeiro, E., and Rossiter, D.: SoilGrids 2.0: producing soil information for the globe with quantified spatial uncertainty, SOIL, 7, 217–240,, 2021. a

Prata, A. J.: A new long-wave formula for estimating downward clear-sky radiation at the surface, Q. J. Roy. Meteorol. Soc., 122, 1127–1151,, 1996. a

Putz, T., Kiese, R., Wollschlager, U., Groh, J., Rupp, H., Zacharias, S., Priesack, E., Gerke, H. H., Gasche, R., Bens, O., Borg, E., Baessler, C., Kaiser, K., Herbrich, M., Munch, J. C., Sommer, M., Vogel, H. J., Vanderborght, J., and Vereecken, H.: TERENO-SOILCan: a lysimeter-network in Germany observing soil processes and plant diversity influenced by climate change, Environ. Earth Sci., 75, 1–14,, 2016. a

Rees, F., Sterckeman, T., and Morel, J. L.: Biochar-assisted phytoextraction of Cd and Zn by Noccaea caerulescens on a contaminated soil: A four-year lysimeter study, Sci. Total Environ., 707, 135654,, 2020. a

Richards, L. A.: Capillary conduction of liquids through porous mediums, Physics, 1, 318–333,, 1931. a, b, c, d

Ritter, A., Hupet, F., Muñoz-Carpena, R., Lambot, S., and Vanclooster, M.: Using inverse methods for estimating soil hydraulic properties from field data as an alternative to direct methods, Agr. Water Manage., 59, 77–96, 2003. a, b

Rummler, T., Wagner, A., Arnault, J., and Kunstmann, H.: Lateral terrestrial water fluxes in the LSM of WRF-Hydro: Benefits of a 2D groundwater representation, Hydrol. Process., 36, e14510,, 2022. a

Schaap, M. G. and Van Genuchten, M. T.: A modified Mualem–van Genuchten formulation for improved description of the hydraulic conductivity near saturation, Vadose Zone J., 5, 27–34, 2006. a

Schaap, M. G., Leij, F. J., and van Genuchten, M. T.: ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions, J. Hydrol., 251, 163–176,, 2001. a, b

Schellekens, J., Dutra, E., Martínez-de la Torre, A., Balsamo, G., van Dijk, A., Sperna Weiland, F., Minvielle, M., Calvet, J.-C., Decharme, B., Eisner, S., Fink, G., Flörke, M., Peßenteiner, S., van Beek, R., Polcher, J., Beck, H., Orth, R., Calton, B., Burke, S., Dorigo, W., and Weedon, G. P.: A global water resources ensemble of hydrological models: the eartH2Observe Tier-1 dataset, Earth Syst. Sci. Data, 9, 389–413,, 2017. a

Schneider, J., Groh, J., Pütz, T., Helmig, R., Rothfuss, Y., Vereecken, H., and Vanderborght, J.: Prediction of soil evaporation measured with weighable lysimeters using the FAO Penman–Monteith method in combination with Richards' equation, Vadose Zone J., 20, e20102,, 2021. a

Schrader, F., Durner, W., Fank, J., Gebler, S., Pütz, T., Hannes, M., and Wollschläger, U.: Estimating Precipitation and Actual Evapotranspiration from Precision Lysimeter Measurements, Proced. Environ. Sci., 19, 543–552,, 2013. a

Schwärzel, K., Šimůnek, J., Stoffregen, H., Wessolek, G., and van Genuchten, M. T.: Estimation of the Unsaturated Hydraulic Conductivity of Peat Soils Laboratory versus Field Data, Vadose Zone J., 5, 628–640,, 2006. a

Séré, G., Ouvrard, S., Magnenet, V., Pey, B., Morel, J. L., and Schwartz, C.: Predictability of the Evolution of the Soil Structure using Water Flow Modeling for a Constructed Technosol, Vadose Zone J., 11, ,, 2012. a, b

Sommer, R. and Stöckle, C.: Correspondence between the Campbell and van Genuchten Soil-Water-Retention Models, J. Irrig. Drain. Eng., 136, 559–562,, 2010. a

Sterckeman, T., Douay, F., Proix, N., and Fourrier, H.: Vertical distribution of Cd, Pb and Zn in soils near smelters in the North of France, Environ. Pollut., 107, 377–389,, 2000. a

Stieglitz, M., Rind, D., Famiglietti, J., and Rosenzweig, C.: An Efficient Approach to Modeling the Topographic Control of Surface Hydrology for Regional and Global Climate Modeling, J. Climate, 10, 118–137,<0118:AEATMT>2.0.CO;2, 1997. a

Tavakoly, A. A., Habets, F., Saleh, F., Yang, Z.-L., Bourgeois, C., and Maidment, D. R.: An integrated framework to model nitrate contaminants with interactions of agriculture, groundwater, and surface water at regional scales: The STICS–EauDyssée coupled models applied over the Seine River Basin, J. Hydrol., 568, 943–958, 2019. a, b

Taylor, K. E.: Summarizing multiple aspects of model performance in a single diagram, J. Geophys. Res.-Atmos., 106, 7183–7192,, 2001. a, b

Tifafi, M., Bouzouidja, R., Leguédois, S., Ouvrard, S., and Séré, G.: How lysimetric monitoring of Technosols can contribute to understand the temporal dynamics of the soil porosity, Geoderma, 296, 60–68,, 2017. a

van Genuchten, M. T.: A Closed-form Equation for Predicting the Hydraulic Conductivity of Unsaturated Soils, Soil Sci. Soc. Am. J., 44, 892–898,, 1980. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Vereecken, H., Maes, J., Feyen, J., and Darius, P.: Estimating the Soil Moisture Retention Characteristic From Texture, Bulk Density, and Carbon Content, Soil Sci., 148, 389–403,, 1989. a, b, c, d

Vereecken, H., Weihermüller, L., Assouline, S., Simunek, Jiri, 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., and Xue, Y.: Infiltration from the Pedon to Global Grid Scales: An Overview and Outlook for Land Surface Modeling, Vadose Zone J., 18, 1–53,, 2019. a, b, c

Vergnes, J.-P. and Decharme, B.: A simple groundwater scheme in the TRIP river routing model: global off-line evaluation against GRACE terrestrial water storage estimates and observed river discharges, Hydrol. Earth Syst. Sci., 16, 3889–3908,, 2012. a

Vergnes, J.-P., Decharme, B., and Habets, F.: Introduction of groundwater capillary rises using subgrid spatial variability of topography into the ISBA land surface model, J. Geophys. Res.-Atmos., 119, 11065–11086, 2014. a

Vergnes, J.-P., Roux, N., Habets, F., Ackerer, P., Amraoui, N., Besson, F., Caballero, Y., Courtois, Q., de Dreuzy, J.-R., Etchevers, P., Gallois, N., Leroux, D. J., Longuevergne, L., Le Moigne, P., Morel, T., Munier, S., Regimbeau, F., Thiéry, D., and Viennot, P.: The AquiFR hydrometeorological modelling platform as a tool for improving groundwater resource monitoring over France: evaluation over a 60-year period, Hydrol. Earth Syst. Sci., 24, 633–654,, 2020.  a, b

Vincendon, B., Édouard, S., Dewaele, H., Ducrocq, V., Lespinas, F., Delrieu, G., and Anquetin, S.: Modeling flash floods in southern France for road management purposes, J. Hydrol., 541, 190–205,, 2016. a

Vogel, H.-J.: Scale Issues in Soil Hydrology, Vadose Zone J., 18, 190001,, 2019. a

Vogel, T., van Genuchten, M., and Cislerova, M.: Effect of the shape of the soil hydraulic functions near saturation on variably-saturated flow predictions, Adv. Water Resour., 24, 133–144,, 2000. a, b, c, d, e

Wafa, M., Boote, K., Cavero, J., and Dechmi, F.: Adapting the CROPGRO model to simulate alfalfa growth and yield, Agron. J., 110, 1777–1790,, 2018. a

Weihermüller, L., Siemens, J., Deurer, M., Knoblauch, S., Rupp, H., Göttlein, A., and Pütz, T.: In Situ Soil Water Extraction: A Review, J. Environ. Qual., 36, 1735–1748,, 2007. a, b

Weynants, M., Vereecken, H., and Javaux, M.: Revisiting Vereecken Pedotransfer Functions: Introducing a Closed-Form Hydraulic Model, Vadose Zone J., 8, 86–95,, 2009. a, b, c, d

Wolf, D. D., Kimbrough, E. L., and Blaser, R. E.: Photosynthetic Efficiency of Alfalfa with Increasing Potassium Nutrition1, Crop Sci., 16, 292–294,, 1976. a

Wösten, J., Pachepsky, Y., and Rawls, W.: Pedotransfer functions: bridging the gap between available basic soil data and missing soil hydraulic characteristics, J. Hydrol., 251, 123–150,, 2001. a, b

Wösten, J. H. M. and van Genuchten, M. T.: Using Texture and Other Soil Properties to Predict the Unsaturated Soil Hydraulic Functions, Soil Sci. Soc. Am. J., 52, 1762–1770,, 1988. a, b, c, d, e, f

Zacharias, S., Bogena, H., Samaniego, L., Mauder, M., Fuß, R., Pütz, T., Frenzel, M., Schwank, M., Baessler, C., Butterbach-Bahl, K., Bens, O., Borg, E., Brauer, A., Dietrich, P., Hajnsek, I., Helle, G., Kiese, R., Kunstmann, H., Klotz, S., Munch, J.C., Papen, H., Priesack, E., Schmid, H. P., Steinbrecher, R., Rosenbaum, U., Teutsch, G., and Vereecken, H.: A network of terrestrial environmental observatories in Germany, Vadose Zone J., 10, 955–973,, 2011. a

Short summary
Seven instrumented lysimeters are used to assess the simulation of the soil water dynamic in one land surface model. Four water potential and hydraulic conductivity closed-form equations, including one mixed form, are evaluated. One form is more relevant for simulating drainage, especially during intense drainage events. The soil profile heterogeneity of one parameter of the closed-form equations is shown to be important.