Impacts of soil management and climate on saturated and near-saturated hydraulic conductivity: analyses of the Open Tension-disk Infiltrometer Meta-database (OTIM)

Saturated and near-saturated soil hydraulic conductivities K_h (mm h⁻¹) determine the partitioning of precipitation into surface runoff and infiltration and are fundamental to soils' susceptibility to preferential flow. Recent studies found indications that climate factors influence K_h, which is highly relevant in the face of climate change. In this study, we investigated relationships between pedoclimatic factors and K_h and also evaluated effects of land use and soil management. To this end, we collated the Open Tension-disk Infiltrometer Meta-database (OTIM), which contains 1297 individual data entries from 172 different publication sources. We analysed a spectrum of saturated and near-saturated hydraulic conductivities at matric potentials between 0 and 100 mm. We found that methodological details like the direction of the wetting sequence or the choice of method for calculating infiltration rates to hydraulic conductivities had a large impact on the results. We therefore restricted ourselves to a subset of 466 of the 1297 data entries with similar methodological approaches. Correlations between K_s and K_h at higher supply tensions decreased especially close to saturation, indicating a different flow mechanism at and very close to saturation than towards the dry end of the investigated tension range. Climate factors were better correlated with topsoil near-saturated hydraulic conductivities at supply tensions ≥ 30 mm than soil texture, bulk density and organic carbon content. We find it most likely that the climate variables are proxies for soil macropore networks created by the respective biological activity, pedogenesis and climate-specific land use and management choices. Due to incomplete documentation in the source publications of OTIM, we were able to investigate only a few land use types and agricultural management practices. Land use, tillage system and soil compaction significantly influenced K_h, with effect sizes appearing comparable to the ones of soil texture and soil organic carbon. The data in OTIM show that experimental bias is present, introduced by the choice of measurement time relative to soil tillage, experimental design or data evaluation procedures. The establishment of best-practice rules for tension-disk infiltrometer measurements would therefore be helpful. Future studies are needed to investigate how climate shapes soil macropore networks and how land use and management can be adapted to improve soil hydraulic properties. Both tasks require large numbers of new measurement data with improved documentation on soil biology and land use and management history.

1 Introduction

Climate models predict more frequent extreme weather events such as high-intensity rainfall with the onset of global warming. To prevent water runoff and erosion, soils need sometimes to be able to conduct large amounts of water in a short time. It is generally accepted that one key soil property is the saturated hydraulic conductivity K_s (mm h⁻¹), as it determines the partitioning of precipitation into surface runoff and infiltration. A large K_s reduces erosion risks and allows water to infiltrate into deeper soil layers, where it may replenish an important reservoir of plant-available water or contribute to groundwater recharge. The hydraulic conductivity of a soil decreases with decreasing water content, i.e. with decreasing water saturation. The hydraulic conductivity in the so-called near-saturated range (between 0 and 100 mm matric tensions) is likewise important. For rainfall intensities smaller than but close to K_s, soils with a larger near-saturated hydraulic conductivity K_h (mm h⁻¹) will remain less water-saturated because they are able to conduct the precipitation water within smaller macropores. Therefore, they are less susceptible to preferential flow (Larsbo et al., 2014) by which agrochemicals and other solutes quickly leach towards the groundwater. Moreover, a large K_h also indicates a well-aerated soil, which drains more quickly and helps air to escape the soil in the case of heavy rainfall. This further reduces the risk of surface runoff and erosion as entrapped air strongly decreases soil hydraulic conductivity.

Saturated hydraulic conductivity is measured either in the laboratory on small cylinders, usually with diameters < 7 cm (Klute and Dirksen, 1986), or is acquired from field measurements using either single- or double-ring infiltrometer methods (Angulo-Jaramillo et al., 2000). In addition, near-saturated hydraulic conductivities can be measured using a tension-disk infiltrometer. The method is designed as a field method but has occasionally been applied in the laboratory. Using a tension-disk infiltrometer, hydraulic conductivities at supply tensions between ca. 0.5 and ca. 60 to 150 mm can be obtained, depending on the specifications of the infiltrometer. All measurement techniques for saturated and near-saturated hydraulic conductivity are laborious, time-consuming and constrained to a relatively small soil volume.

It is necessary to develop pedotransfer functions to estimate soil hydraulic conductivities for large-scale modelling applications, as we cannot measure everywhere (Bouma, 1989; Van Looy et al., 2017; Wösten et al., 2001). The development of a pedotransfer function requires a database from which it can be derived. For example, the well-known pedotransfer function ROSETTA (Schaap et al., 2001) is based on the open UNSODA database (Nemes et al., 2001). The equations published in Tóth et al. (2015) are derived from the proprietary EU-HYDI database (Weynants et al., 2013). The pedotransfer functions of Jarvis et al. (2013) are based on an unpublished meta-database containing tension-disk infiltrometer data. Collecting published measurements of saturated and near-saturated hydraulic conductivity measurements into meta-databases and pairing them with other existing databases is essential for developing pedotransfer functions. A notable example is the SWIG database (Rahmati et al., 2018) that collates more than 5000 data sets from soil infiltration measurements covering the entire globe. Another big effort in collecting information on saturated hydraulic conductivity is the newly published SoilKsatDB (Gupta et al., 2021a), which combines saturated hydraulic conductivity data from several large databases, amongst others UNSODA and SWIG, with additional measurements published in independent scientific studies. However, none of the databases cited above provides open-access infiltration measurements at tensions near saturation (h > 0 mm), which limits their use to the estimation of saturated hydraulic conductivity.

While reasonably good estimations of K_s from easy-to-measure or readily available site properties appear to be possible for peat soils (Morris et al., 2022), pedotransfer functions for K_s of mineral soils exhibit poor predictive performance, with coefficients of determination R² not exceeding 0.25 (Weynants et al., 2009; Jorda et al., 2015). Early approaches, like HYPRES (Wösten et al., 1999) and ROSETTA, focused solely on soil properties like texture, bulk density and organic carbon content as predictors of K_s. At the time, it was not sufficiently recognized that soil K_s is mostly determined by the morphology of macropore networks, especially in finer-textured soils (Vereecken et al., 2010; Koestel et al., 2018; Schlüter et al., 2020). A pedotransfer function for K_s therefore ideally requires a database that contains direct information on the macropore network itself. However, since such measures are even more cumbersome and time-consuming to obtain (e.g. by X-ray tomography) than measuring hydraulic conductivity itself, it is more reasonable and makes more sense to use proxies from which the macrostructure in a soil can be inferred. Ideal candidates would be root growth and the activity of soil macrofauna, which both strongly determine the development of macropore networks in soil (Meurer et al., 2020). However, they are difficult to measure. Proxies that are more promising are land use and farming practices, such as tillage or soil compaction due to trafficking. Plant growth and soil macrofauna in turn are influenced by the local climate. The climate also sets boundaries for the land use and the associated soil management practices and thus provides feedback on root growth and macrofaunal activity. Wetting and drying cycles and thus the formation and closure of cracks are also regulated by the climate, as are splash erosion and soil crusting. It is therefore not surprising that climate variables typically are correlated with saturated and near-saturated hydraulic conductivities (Jarvis et al., 2013; Jorda et al., 2015; Hirmas et al., 2018; Gupta et al., 2021b). Jorda et al. (2015) found that land use itself was the most important predictor of saturated hydraulic conductivity.

The time of measurement of the hydraulic conductivity (or soil sampling) also has a crucial impact. In an agricultural soil, the hydraulic properties of a freshly prepared seedbed differ from those measured later at harvest. Several studies have demonstrated the evolution of hydraulic conductivity with time (Messing and Jarvis, 1990, 1993; Bodner et al., 2013; Sandin et al., 2017). Soil management options (such as tillage or the use of cover crops) actively influence the soil-saturated and near-saturated hydraulic conductivity. Information on their impact is therefore especially important but so far has hardly been investigated in meta-studies.

In this study, we focused on quantifying the effect of soil management practices on soil-saturated and near-saturated hydraulic conductivity, K_h. We also investigated relationships between K_h and other potentially important influencing factors like soil properties, local climate and details of the measurement methods. In this process, we expanded and published the previously unpublished meta-database of tension-disk infiltration measurements that was first reported by Jarvis et al. (2013). We refer to this database as OTIM in the following (Open Tension-disk Infiltrometer Meta-Database). It complements the currently available public databases on hydraulic conductivities, which are strongly based on laboratory measurements or ring infiltrometer methods.

2 Material and methods

2.1 Meta-database, OTIM

2.1.1 Data collection

The first version of OTIM was compiled for the study by Jarvis et al. (2013). The original database contained 753 tension-disk infiltrometer data entries collated from 124 source publications covering 144 different locations around the globe. We have extended this database by 544 new tension-disk infiltrometer data entries from 48 additional studies that were published after 2012. The search for publications was carried out between 31 May 2021 and 23 June 2021 using the queries and search engines detailed in Table A1.

Table 1 List of new entries added to the Jarvis et al. (2013) database.

We found 115 publications containing tension-disk infiltrometer measurements published in 2013 or later. We retained the data for further analysis when (i) K_h or the infiltration rate was measured at more than two tensions larger or equal to 5 mm and (ii) sufficient meta-data on soil and site properties (at least soil texture) and soil management practices (at least land use and tillage) were available. If a publication only reported infiltration rates, we calculated hydraulic conductivity using the method of Ankeny et al. (1991). Only 45 of the 115 publications fulfilled the above-mentioned criteria. Table A2 summarizes how many papers were rejected and for which reasons. For 27 of the 45 retained studies, we digitized the published K_h values from figures using WebPlotDigitizer (open-source web-based software created by Ankit Rohatgi, https://automeris.io/WebPlotDigitizer/, last access: 14 July 2023). For cases in which K_h measurements were mentioned in a publication but the results were not reported, we contacted the corresponding authors. We received the data in this fashion for three of these publications (Alletto et al., 2015; Larsbo et al., 2016a; Meshgi and Chui, 2014). The new studies containing data were added to the OTIM database and are summarized in Table 1.

In addition to adding data from new publications to OTIM, we revisited the studies contained in the original version of the database and collected additional information on soil management practices associated with the measured data. For each soil management option, OTIM contains two columns. In the first column, the information as given in the source publication is stored. The second column summarizes this information into a few classes, which were subsequently used in the meta-analysis. In this study, we investigated effects of land use, tillage system, soil compaction and day of measurement relative to the latest tillage operation on the field. A compaction class was assigned to a data entry only if the plot had been described as “compacted” or “not compacted” in the source publication. Compacted data entries corresponded, for example, to infiltration measurements in wheel tracks or to plots of a compaction experiment. The day of measurement relative to tillage was also included, with the data labelled “freshly tilled” when the authors in the source publication stated that the measurements had taken place soon after tillage. Otherwise, it was assumed that the soil already had time to consolidate before the infiltration measurements were carried out. All soil texture data were mapped onto the USDA classification system using the method proposed in Nemes et al. (2001).

2.1.2 Climate data and soil classification

The climatic data entries provided in the database were created using the bioclimatic raster data (BioClim) provided by WorldClim (https://worldclim.org, last access: 14 July 2023). The data were averaged across the years 1970 to 2000 and had a 30 arcsec resolution (∼ 1 km²; Fick and Hijmans, 2017). The available climate variables were mean annual temperature, mean annual precipitation, the mean temperature and mean precipitation of the warmest, coldest, wettest and driest quarters and months, respectively, the isothermality, the mean diurnal and annual temperature range, the seasonality for temperature and precipitation. In addition to the bioclimatic data in WorldClim, we included the aridity index (here defined as the annual precipitation divided by the potential evapotranspiration) and the average annual potential evapotranspiration (ET₀). Both were inferred from the Global Aridity Index and Potential Evapotranspiration Climate Database v2 that is based on the WorldClim database (Trabucco and Zomer, 2019). The World Reference Base (WRB) soil type was also extracted from the source publications. When this was not reported, the SoilGrids database by ISRIC (Poggio et al., 2021) was used to infer it. The map contained information about the main soil type regarding the WRB classes (IUSS Working Group WRB, 2015). The most probable soil type was chosen for each location. For all the discussed climate and soil maps, the Python package rasterio (v1.2.10) was used to collect the variables from the corresponding raster cell at the location coordinates given in the source publications.

2.1.3 Model fit to infer K_h at unmeasured near-saturated tensions

Tension-disk infiltrometers measure infiltration rates at a specific supply tension (Angulo-Jaramillo et al., 2000). They consist of a ceramic disk to which a water reservoir and a bubbling tower are attached. The ceramic disk is saturated and hydraulically connected to the soil by inserting a layer of fine sand between the disk and the soil surface. The supply tension at the bottom of the ceramic disk is adjusted by the bubbling tower. The measured unconfined (i.e. three-dimensional) infiltration rates are then commonly converted to hydraulic conductivities with the aid of the Wooding equation (Wooding, 1968). Note that unconfined tension-disk infiltrometers cannot provide measurements at a tension of zero, i.e. K_s. Even if many publications report K_s values obtained from tension-disk infiltrometers, these measurements must have been conducted at tensions slightly larger than zero, as water would otherwise have freely leaked out of the tension disk. For this reason, we set the tensions for K_s measurements to 1 mm but still refer to these data as saturated hydraulic conductivity. Note that we discuss matric potentials in terms of tensions (negative pressures) throughout this paper. For convenience, we denote K_h at a specific tension by replacing the subscript “h” with the tension value in millimetres. For example, K₁₀₀ denotes K_h at a supply tension of 100 mm.

Following Jarvis et al. (2013), we interpolated K_h for tensions between the ones measured in the source publications. We achieved this by fitting a log–log linear model with a kink at a tension h_min, which denotes the tension at which the largest effective pores in the soil are water-filled (see Fig. 1). Therefore, K_h≡K_s for all tensions h≤h_min. If K_s was not measured but instead a K_h value at h≤ 5 mm was available, K_s was set to the available K_h value (Fig. 1, orange line). In cases where more than one K_h value was measured at a tension smaller than or equal to 5 mm (including h = 0 mm, i.e. K_s), we averaged them and fixed K_s and K_h for h≤h_min to the average (Fig. 1, green line). K_h values at h > 5 mm were used to fit the log–log linear relationship. The tension at which the fitted log–log slope intersected with K_s is defined as h_min. We used the fitted model to estimate all K_h values for tensions for $\mathrm{10}\le h\le$ 100 mm at 10 mm intervals. The K_h values were only interpolated between the tensions that were measured in the source publication. The only exceptions to this rule were made in the case where a K_h value for a tension of 80 or 90 mm was provided together with at least one other K_h value measured at a smaller tension. Then, the missing K_h values were extrapolated up to a tension of 100 mm. Figure 1 shows examples of model fits. Only entries with an R² greater than or equal to 0.9 were retained in the analysis.

Figure 1 Two examples of the linear fit in log–log space. The colours denote two differently measured tension series. The filled circles correspond to measured K_h, while the lines indicate the interpolation carried out by the model. The bold black dashed line marks a supply tension of 5 mm. K_h values at tensions between 0 and 5 mm were assumed to be identical to K_s. Reported K_s values were assigned a tension of 1 mm for illustration purposes. The equations for the linear part of the fit are shown in the legend. C represents the intercept with the y axis of the linear fit in log space and h_min corresponds to the supply tension at which the largest pores in the soil are water-filled.

Figure 2 Map of the study locations collected in OTIM. The values are shown for the filtered entries (“focus”) and in parentheses for all the entries available in the database (“focus” and “extra”).

2.2 Data availability and spatial coverage

Although 92 % of the OTIM data are from topsoils, OTIM also contains some data points measured at greater soil depths. In the following meta-analysis, only measurements from the topsoil were included to prevent bias, and all data sets measured at soil depths below 200 mm were removed. Last but not least, we found that the relationship between supply tension and K_h was distorted if data entries were included that did not cover the complete tension range from h = 0 to 100 mm. Possible reasons for the difficulties in matching K_h data from tension series with different lengths are discussed at the beginning of the Results and Discussion sections. Otherwise, we focused on data entries that included K_h values for the complete tension range in the exploratory data analysis and the meta-analyses. The available data sets after these filtering steps correspond to the ones indicated in blue (and termed “focus”) in the following figures.

Most tension-disk infiltrometer studies were conducted in Europe, North America and south-eastern Australia (Fig. 2). Clearly, fewer studies have been carried out in Asia, South America and Africa. The lack of data sets from Russia, Mesoamerica, the Arctic regions and the tropics is remarkable. This geographical bias is aggravated if only measurements of the topsoil are considered that allow inferences of K_h for the complete range of tensions ($\mathrm{0}\le h\le$ 100 mm) with a sufficiently good coefficient of determination. Then, almost all the data entries collected in southern South America and south-eastern Australia were omitted as well. Overall, the data in OTIM mostly stem from temperate climate regions.

Figure 3 (a) Number of available K_h values per supply tension, (b) available tension series with the black bar indicating the span between T_min and T_max and (c) their respective frequency in the database. The values are shown for the filtered entries (focus) and for all the entries available in the database (focus and extra).

Table 2 Number of entries and gaps for each feature along with units and ranges (if continuous) or choices (if categorical). The values are shown for the filtered entries (“focus”) and in parentheses for all the entries available in the database (“focus” and “extra”). CV stands for coefficient of variation.

Figure 3 depicts the number of K_h values available for $\mathrm{0}\le h\le$ 100 mm. These figures represent the hydraulic conductivities derived from the log–log linear model presented above, not the raw data measured and reported in the source publications. A large number of entries spans the full range of tensions of interest (0 to 100 mm), whereas a smaller number of entries only has data up to a tension of 60 mm. Often, but not always, such data series were obtained with the widely available mini-disk infiltrometer distributed by the Meter group (formerly by Decagon), which is limited to tensions h≤ 70 mm. An overview of the meta-data included in OTIM is given in Table 2. Data gaps are present, especially for bulk density and for information on the soil management at the study site, apart from tillage operations. Note that the annual mean temperature and precipitation are only two examples representing the climatic variables enumerated in Sect. 2.3. There are very few missing values for the climate data, since they were estimated from the coordinates of the study sites. The same holds for the elevation data and information on the WRB soil type.

Figure 4 Distributions of continuous and categorical variables in the focus data set.

The meta-data for the data sets used in the exploratory data analysis are summarized in Fig. 4. OTIM contains predominantly data from arable fields. The distributions of the climate variables confirm that the data in OTIM were also mostly acquired in temperate climates, with a bias towards the somewhat drier climates that are most typical for arable land. The soil texture, bulk density and organic carbon content data also appear to be reasonably representative of soils in this climate zone.

2.3 Exploratory data analysis

Some source publications only provided a few data entries for K_h, sometimes only comparing two different treatments, while other source studies contain data for a larger number of treatments and/or sites. In some publications, data for all individual tension-disk measurements are available even if replicates were measured. In others, only averages of the replicated measurements are reported, while still others yield average K_h values for individually replicated treatment blocks. This makes appropriate data weighting complicated but also extremely important when analysing the meta-data set. It also introduces uncertainty, because it is not always clear whether the replicated averages were calculated using the geometric or arithmetic mean. Considering that hydraulic conductivities at or near saturation are known to be log-normally distributed, the former would be best. In the following, we assumed that geometric averaging was used when replicated values were reported in source publications. In the following, we calculated data weights as

$$\begin{array}{}\text{(1)}& {\mathit{\omega }}_i={\displaystyle \frac{n_{\mathrm{r},i}}{\sqrt{N_i}}},\end{array}$$

where ω_i is the weight for data entry i, n_r,i is the number of replicates from which the values of i were averaged, and N_i is the total number of measurements included in the publication from which data entry i was obtained. With this approach, we up-weighted data entries according to the number of replicate measurements from which they were averaged and down-weighted the impact of studies that published larger numbers of data.

We used weighted Spearman rank correlation coefficients to investigate relationships between continuous variables. We considered correlations to be significant if they exhibited p values of less than 0.05. The latter were determined numerically by running randomization tests with 200 repetitions.

Table 3 Number of studies and paired comparisons with their respective control and treatment values used for the meta-analysis exemplary for the K₁₀₀ values.

2.4 Meta-analysis

Data entries in OTIM with specific land use or management were very unevenly distributed. For example, the large majority of data were measured on sites with land use “arable” (see Fig. 4a). Such uneven distributions may lead to bias when averaged over all entries of a specific feature in exploratory data analyses. We therefore investigated the effects of land use and management as well as soil compaction and time of measurement on K_h with the aid of pairwise comparisons published within individual studies and calculated effect sizes (ESs) for each investigated class.

To reduce bias arising from the varying number of data entries published within individual studies, we grouped all entries according to the factors land use, tillage, compaction and sampling time. Here we only considered binary pairs, which are arable or not arable in the case of land use and tilled or not tilled, compacted or not compacted and “measured soon after tillage” or “measured on consolidated soil” for the other three factors. In addition, we checked whether different entries within individual studies stemmed from the same site or a very similar site. We did this by comparing the respective USDA texture classes and a climate variable, i.e. the aridity class. All data entries within each individual study that exhibited identical land use, soil management, soil compaction, sampling time, texture and aridity were averaged, and the number of corresponding replicates was summed.

For each binarized factor (e.g. tillage), a control value was chosen (e.g. zero tillage). All values different from the control were represented by the treatment (e.g. conventional tillage and reduced tillage). Within individual studies, pairs among the averaged entries were formed for each combination of a control value and a treatment value. These pairs were used to compute the effect size. Following Basche and DeLonge (2019), we defined the effect sizes as log₁₀ of the ratio of K_h,t of the treatment divided by K_h,c of the control:

$$\begin{array}{}\text{(2)}& {\mathrm{ES}}_l={\text{log}}_{\mathrm{10}}\left({\displaystyle \frac{K_{\mathrm{h},\mathrm{t}}}{K_{\mathrm{h},\mathrm{c}}}}\right),\end{array}$$

where the subscript l indicates the lth pair for which the effect size was computed and the indices “t” and “c” stand for the treatment and control, respectively. The average ES for each of the four investigated factors was calculated as the weighted mean of the individual ES_l using the weight

$$\begin{array}{}\text{(3)}& w_l={\displaystyle \frac{{\mathit{\nu }}_{\mathrm{c}}{\mathit{\nu }}_{\mathrm{t}}}{{\mathit{\nu }}_{\mathrm{c}}+{\mathit{\nu }}_{\mathrm{t}}}},\end{array}$$

where the subscript l indicates again the lth pair for which the effect size was computed and ν_c and ν_t denote the number of (summed) replicates for the control and treatment, respectively. In addition, we calculated the weighted standard error

$$\begin{array}{}\text{(4)}& {\mathit{\sigma }}_{\stackrel{\mathrm{‾}}{\mathrm{ES}}}=\sqrt{{\displaystyle \frac{{\sum }_{i=\mathrm{0}}^n{w}_l{\left({\mathrm{ES}}_l-\stackrel{\mathrm{‾}}{\mathrm{ES}}\right)}^{\mathrm{2}}}{\frac{n-\mathrm{1}}{n}{\sum }_{l=\mathrm{0}}^n{w}_l}}},\end{array}$$

where $\stackrel{\mathrm{‾}}{\mathrm{ES}}$ is the mean effect size. Table 3 summarizes the evaluated factors, the number of pairs involved and the number of different studies from which the pairs were obtained.

To estimate the robustness of the effect size, we carried out a sensitivity analysis using the jackknife technique, similarly to Basche and DeLonge (2019). This method aims to show the sensitivity of the averaged effect size to data from specific studies. For each factor, a given number of studies was randomly picked and removed from the data set. The averaged effect size and its standard error were computed with the rest of the data set. The process started by removing one study, after which up to nine more studies were removed. This random selection was repeated 50 times to rule out bias. The average of the means and standard errors for the 50 realizations was computed and plotted. Observed effect sizes were judged to be trustworthy if they did not change after removal of studies to calculate them. We constrained the sensitivity analyses in our study to the effect sizes for K_s and K₁₀₀.

Figure 5 Evolution of the weighted mean K_h with tension available in OTIM, sorted by the tension range the data were spanning. The number of publications from which the data originated is shown in parentheses in the legend. The shaded areas and the error bars represent the weighted standard error of the mean.

Figure 6 Evolution of the weighted mean K_h as a function of applied tension for the (a) disk diameter, (b) direction and (c) method of fitting. “S.-S.” stands for “steady-state”. More specifically, the method “S.-S. constant” is outlined in Logsdon and Jaynes (1993), “`S.-S. multi-disc” in Smettem and Clothier (1989), “S.-S. piece-wise” in Reynolds and Elrick (1991) or Ankeny et al. (1991) and “Transient” in Zhang (1997) or Vandervaere et al. (2000). The shaded areas and the error bars represent the weighted standard error of the mean. The bar plots in each subplot indicate how many data points of each class were in the data set.

3 Results and discussion

3.1 Differences between data entries with different tension ranges

If all data are considered (focus and extra), Fig. 3 illustrates that approximately 40 % of the data in OTIM provided K_h for every h with $\mathrm{0}\le h\le$ 100 mm. For another 40 %, K_h was only measured in the wet range, i.e. at tensions below 70 mm. The remaining K_h data were only acquired in the dry range. Here, we counted all data entries for which K_s was not measured and could not be estimated. Figure 5 shows how data from entries with complete, dry and wet ranges differed. The K_h for the wet range receded more quickly with increasing tension than series that also included measurement in the dry range. A large portion of these data sets was obtained with the mini-disk infiltrometer. However, a closer inspection of the impact of the disk diameters used to acquire the respective K_h did not confirm suspicions that the bias was related to the use of this special type of infiltrometer (see Fig. 6a). Instead, the observed differences between the K_h curves could have been introduced by co-correlations with soil texture or climate. Another explanation may be experimenter bias, since individual research groups tend to use specific tension ranges for more than one study. In this study, however, we focused solely on data entries for which we were able to reconstruct K_h for all h between 0 and 100 mm supply tension in the following exploratory data analyses and meta-analyses. This greatly facilitated the data interpretation.

Figure 7 Weighted Spearman rank correlation coefficients between K_h at different tensions. Correlation coefficients are shown up to p values of 0.001.

3.2 Statistical relationships between K_h and the methods used

Figure 6a confirms that the diameter of the tension disk did not have a systematic impact on the results. The majority of the data were collected starting under dry conditions (large tensions) and subsequently measured under increasingly wet conditions (smaller tensions). Figure 6b illustrates that beginning the experiment under wet conditions is associated with larger hydraulic conductivities at identical supply tensions. This is well known and is referred to as hysteresis, which is due to ink-bottle effects, impacts of water repellency, air entrapment and swelling of clay particles (Hillel, 2003). Figure 6c shows that the large majority of studies used the “steady-state piecewise” method to solve the Wooding equation and convert the measured infiltration rates to hydraulic conductivities. This method leads to smaller K_h for larger tensions than the other methods. The “transient” and “steady-state constant” methods yielded a larger K_h in the unsaturated range. For the latter method, it is known that it overestimates unsaturated K_h (Jarvis et al., 2013). We tested whether excluding data from transient and steady-state-constant methods changed the results of the meta-analyses but found that they only changed to a minor degree. Data from all the methods were therefore included in the following. Note that the transient method was mostly applied in conjunction with mini-disk infiltrometers, but the respective data are not included in Fig. 6 since they do not span the entire suction range.

3.3 Correlation between K_h at different tensions

The fact that correlations between K_h estimated at supply tensions between 40 and 100 mm were relatively stable (Fig. 7) indicates that the respective flow paths and/or mechanisms remained very similar in this tension range. However, these correlations weakened at tensions between 10 and 20 mm, giving rise to the existence of a threshold above which water flow in the largest macropores becomes the dominant flow mechanism. The poor correlation between K_s and K_h at larger supply tensions is in line with findings that K_s is not well suited for inferring soil-unsaturated hydraulic conductivities (Schaap and Leij, 2000).

Figure 8 Evolution of the weighted mean K_h as a function of applied tension for (a) soil texture, (b) soil bulk density and (c) soil organic carbon. The shaded areas and the error bars represent the weighted standard error of the mean. The soil textures were classified using USDA texture classes as follows: fine (clay, clay loam, silty clay, silty clay loam), medium (silt loam, loam) and coarse (loamy sand, sand, sandy clay, sandy clay loam, sandy loam). The bar plots in each subplot indicate how many data points of each class were in the data set.

3.4 Statistical relationships between K_h and soil properties

Soils with coarse texture exhibited larger K_h in the unsaturated range, which is caused by the large and abundant primary pores in between individual sand grains (Fig. 8a). At saturation, the average hydraulic conductivity of all three texture classes was similar. This is explained by the presence of large structural pores in the medium- and fine-textured soils. Medium-textured soils had the lowest K_h in the investigated range of tensions, which may be due to a denser soil matrix in loamy soils and a lower structural stability of silty soils. Larger bulk densities decreased K_h across the whole range of investigated tensions, which reflects the reduced porosity with increasing bulk density (Fig. 8b).

Figure 9 Weighted Spearman rank correlation coefficients between climate variables, elevation above sea level and soil properties. Correlation coefficients are shown up to p values of 0.001.

The hydraulic conductivity in the saturated and near-saturated ranges is especially affected by soil compaction, which predominantly reduces the abundance and connectivity of macropores (Pagliai et al., 2004; Whalley et al., 1995). Large bulk densities are also known to reduce burrowing activities of the soil macrofauna (Capowiez et al., 2021) and root growth (Lipiec and Hatano, 2003), also leading to less abundant and less connected large macropores. An increase in the soil organic carbon (SOC) content was connected to smaller K_h at the dry end of the investigated tension range if soils with organic carbon contents of more than 0.03 kg kg⁻¹ were excluded (Fig. 8c). This decrease may be explained by water repellency, which is generally positively correlated with organic carbon content. A similar observation was already reported in Jarvis et al. (2013). Note that no major correlations of SOC with soil texture were observed in the investigated data set (Fig. 9). For soils with organic carbon contents larger than 0.03 kg kg⁻¹, K_h increased once again. This may indicate that, above this threshold, better-developed macropore networks associated with large SOC contents (e.g. Larsbo et al., 2016b) outweighed any effects of water repellency.

Figure 10 Weighted Spearman rank correlation coefficients between K_h at different tensions, climatic features and soil properties. Correlation coefficients are shown up to p values of 0.001.

3.5 Statistical relationships between K_h and climate variables

One important observation made in recent years was that saturated and near-saturated hydraulic conductivities correlated strongly with climate variables (Jarvis et al., 2013; Jorda et al., 2015; Hirmas et al., 2018). Figure 10 gives an overview of weighted Spearman rank correlations between K_h and 6 of the 20 climate variables included in OTIM that exhibited the strongest correlations with K_h. The elevations of the sampling site above sea level together with its latitude and longitude, soil texture, bulk density and soil organic carbon content are also shown for comparison. It is striking that the soil properties were less well correlated with K_h than some of the climate variables. Of the three USDA texture fractions, the clay content was negatively correlated, the sand content was positively correlated with K_h in the drier investigated tension range, and no significant correlations were found for the silt fraction (Fig. 10). Only the bulk density exhibited correlation coefficients as large as the climate variables.

The largest absolute values of the weighted rank correlations were observed for the mean diurnal range of temperature and the aridity index. Both reach a maximum at the dry end of the considered tension range, i.e. for K₁₀₀, with correlation coefficients of 0.43 and −0.4, respectively. Figure 9 reveals that both of these best-correlated climate variables were accidently correlated with choices in experimental design and data evaluation made by the investigators in the respective source studies, which will amplify these observed correlations with K_h. However, if a smaller data set is considered in which such methodological bias and potential bias due to differences in land use were eliminated, the correlations persist (Fig. B1). We therefore infer that the observed effect of climate on K_h is real.

The annual mean diurnal temperature range and the aridity index were strongly correlated with each other, with a weighted correlation coefficient of −0.68 (Fig. 9). Strong correlations with at least one of these two variables with absolute values > 0.6 were also found for most of the investigated climate variables. It is difficult to separate the climate effects due to these strong inter-correlations. Nevertheless, it is striking that the mean annual diurnal temperature ranges are much better correlated with K₁₀₀ than the mean annual temperature itself (Fig. 10). In addition, the mean annual precipitation in the driest quarter of the year and the precipitation in the driest quarter of the year exhibited stronger correlations than the mean annual precipitation. It appears that temperature and precipitation fluctuations are more strongly coupled to near-saturated hydraulic conductivities than the absolute temperatures or precipitation amounts.

Among possible reasons for the observed correlations may be increased splash erosion during heavy rainfalls that are common in regions with large precipitation seasonality, more soil compaction in wetter climates due to trafficking, a larger vertical burrowing activity of soil fauna in climates with large diurnal temperature ranges, more vertically oriented root systems in arid climates, or climate-specific choices in land use and soil management. The data in OTIM cannot provide an answer to these questions. Investigations of such relationships should be the focus of future studies.

Another site factor that is positively correlated with K_h is the elevation above sea level (Fig. 10). Notably, elevation above sea level was also found to be an important predictor for K_s in Gupta et al. (2021b), which suggests that there are indeed pedogenetic reasons behind the observed correlation. In the case of infiltrometer measurements, the decreased atmospheric pressure with height on the supply tension can be neglected. The supply tension is always equivalent to the weight of the water column adjusted in the bubbling tower. The weight of the water column will be smaller due to the general decrease in earth's gravitational constant with height due to a larger centrifugal force. However, the weight of the water column would only be reduced by approximately 1 or 2 %. Also, indirect influences of larger heights on the infiltration rate cannot explain the observed correlation. A lower temperature would make the water column denser. However, the effect would be less than 1 % in the relevant temperature range. In contrast, a lower temperature would increase the viscosity of water to a much larger degree, e.g. by up to approximately 30 % between temperatures of 10 and 20 ^∘C. The temperature effect should thus lead to a negative correlation between elevation and K_h, which is the opposite of what was observed. Bias in the K_h measurements due to such physical effects can thus be ruled out. Elevation may instead be a proxy for well-drained soils, as stagnant soil water and high groundwater tables are less likely with height above sea level. This may favour soil life and better-developed root systems and decrease risks of compaction when the soil is trafficked.

The observed correlations of K_h with latitude and longitude probably reflect co-correlations with climate variables together with experimenter bias, since it appears likely that approaches in setting up tension-disk infiltrometers systematically vary between continents, e.g. America and Europe.

Bulk density was the only soil property that exhibited a (negative) correlation strength of > 0.3 with any K_h (Fig. 10). The underlying reasons were discussed above. Notably, the strongest correlations were found at and very close to saturation, probably due to the detrimental effect of soil compaction on macroporosity and macropore connectivity. The more compaction, the less macroporosity and the higher the bulk density, which in turn decreases root growth and bioturbation (Capowiez et al., 2021; Lipiec and Hatano, 2003). The only pedoclimatic factor with a relatively large correlation strength (−0.24) with the saturated hydraulic conductivity was the bulk density (Fig. 10).

Figure 11 Weighted mean log₁₀ response ratio (effect size) of K_h from K₁₀₀ to K_s for different management practices where the controls were “not arable”, “no tillage”, “no compaction” and “consolidated soil”, respectively. A positive effect size means that the value of the treatment is greater than the control. The dashed line shows “no effect” (no difference between the treatment and the control). Error bars represent the weighted standard error of the mean.

Figure 12 Sensitivity analysis of the weighted effect size of K at 100 mm tension and K_s for the management practice investigated using the jackknife technique. The error bars represent the standard error.

3.6 Effects of land use, tillage, compaction and sampling time

The average log₁₀ response ratios shown in Fig. 11 illustrate the effects of land use and soil management on K_h for $\mathrm{0}\le h\le$ 100 mm. Note that a value of ± 0.3 in the log₁₀ response ratio corresponds to a factor of 2. Hence, K_s for uncompacted soil was found to be approximately twice as large as for compacted soil (see Fig. 11c). Arable land exhibited a clearly smaller K_s than grasslands and forests, which is in line with observations made by Basche and DeLonge (2019). This difference became smaller with higher tensions (Fig. 11a). The large difference in K_h close to saturation was likely related to traffic compaction and tillage operations that were applied to the majority of the investigated arable soils, which led to the destruction of connected biopores and hence a reduced K_s. On the other hand, tillage breaks up intact soil into individual soil aggregates, which creates, at least initially, a well-connected network of inter-aggregate pores that increase K_h in the near-saturated range (Sandin et al., 2017; Schlüter et al., 2020). This effect of tillage can explain why near-saturated K_h under conventional and reduced tillage was larger than under no tillage (Fig. 11b). However, in this case, even K_s was larger in the tilled fields. It is likely that K_s was reduced in the no-tillage treatments due to traffic compaction on the fields and a lack of soil loosening by tillage as compared to conventionally tilled treatments. Note however that we only investigated topsoils in this study. It is not clear how different tillage types affect K_h in the subsoil. The impact of soil compaction on K_h was clearly negative in the entire investigated range of tensions (Fig. 11c), which is explained by the reduction in porosity and especially the macroporosity during compaction (see also Fig. 8b). In contrast, if the K_h measurements were carried out shortly after tillage operations, K_h was increased for all investigated tensions, especially very close to saturation (Fig. 11d). This confirms that tillage initially increases K_s but that subsequent soil consolidation preferentially disconnects the largest macropores. As a consequence, K_h at and very close to saturation is reduced more strongly than K_h for higher tensions (see Fig. 11a).

Figure 13 Evolution of the weighted mean K_h as a function of applied tension for the (a) aridity class, (b) potential annual evapotranspiration and (c) mean annual diurnal temperature range. The shaded areas and the error bars represent the weighted standard error of the mean.

Figure 12 shows the results of the sensitivity analyses for the effect sizes depicted in Fig. 11. The effect of land use for K₁₀₀ turned out to be most sensitive to the removal of studies (Fig. 12a). The direction of the effect even changed after the removal of six studies, indicating that higher K₁₀₀ values for arable compared to non-arable fields were not just occasional observations but occurred more frequently. More studies would be needed to properly characterize the effect of land use on K₁₀₀. The remaining sensitivity analyses for all the other factors showed that removal of studies did not change or destabilize the results for both K₁₀₀ and K_s.

3.7 Comparison of the effect size of land use, management and sampling time with the effect of climate and soil properties

Effect sizes could only be computed for land use and management, compaction and sampling time. It is therefore difficult to relate the impact of these factors to the ones of measurement method, climate variables and soil properties. Comparisons between Figs. 6, 8 and 13, on the one hand, with Fig. 11 provide some insight. Land-use- and management-related effects together with sampling time (Fig. 11) seem to have a similar effect on K_h to the soil properties (Fig. 8) and the measurement method (Fig. 6). Climate variables seem to have a larger impact on K_h at the dry end of the investigated tension range but a smaller one close to saturation (Fig. 13).

4 Conclusions

Our results suggest that climate change will influence soil hydraulic properties near saturation. This may complicate model predictions of water balance in a future climate, particularly the risks of surface runoff, soil erosion and waterlogging. Climatic factors are more strongly correlated with near-saturated hydraulic conductivities than soil texture, bulk density and organic carbon content. At and very close to soil saturation, the correlations between hydraulic conductivity and climate variables vanished. Instead, the soil bulk density showed the largest correlation, in line with the fact that more compact soils tend to lack a well-connected macropore system. Hypotheses as to why climate variables are correlated with the hydraulic conductivity were discussed but need to be investigated in future studies. Most probably, the impacts of climate are linked to macropore networks associated with biological activity, pedogenesis and land use. Only a few land-use- and soil-management-related factors could be investigated in our study. They were all found to significantly influence K_h, with an effect on sizes similar to those of soil properties like texture and organic carbon content. Also, experimenter bias as introduced by the choice of measurement time relative to soil tillage, experimental design or data evaluation appeared to be as important for the saturated and near-saturated hydraulic conductivities as soil texture or bulk density. There is a need for better documentation and accessibility of measurement data and the associated meta-data, as has already been suggested by others (McBratney et al., 2011; Basche and DeLonge, 2019). OTIM offers the possibility of deriving more comprehensive pedotransfer approaches than the ones in Jorda et al. (2015). The construction and evaluation of such pedotransfer functions are envisioned for an upcoming companion paper to this study.

Appendix A

A1 Data query details

Table A1 Query strings, search engines, number of result pages that were processed and dates of the search for new data for OTIM.

A2 Data rejection

Table A2 Reasons for data rejection.

A3 Database organization

OTIM is organized into nine individual tables illustrated in Fig. A1. The main table is named experiments. It contains identifiers with which all the other tables are linked. The identifiers are shown in bold font in Fig. A1. The reference table contains information on the references for each study. The location table lists the coordinates of the measurement sites. The tables soilProperties, soilManagement and climate store data, as implied by their names. The method table gives details of the specifications of the tension-disk infiltrometer and the method to calculate hydraulic conductivity from the infiltration rate. The rawData table contains the hydraulic conductivities and respective supply tensions as stated in the corresponding source publication. Note that OTIM does not contain raw data for the entries of the original version compiled for Jarvis et al. (2013). Finally, modelFit reports K_h for $\mathrm{0}\le h\le$ 100 mm as described above. For more details, the reader is directed to the “Description” tab of the database (not shown in Fig. A1), where the meanings and units of each column are explained.

Figure A1 Structure of the OTIM database with its different tables and columns. In the soilManagement table, the columns with the suffix “Class” denote columns in which the data reported in the source publications were summarized into classes to facilitate comparisons between them. For example, the reported CurrentCrop like wheat, rye, barley or oat was assigned the CropClass cereals. The rows in bold denote unique identifiers with which the table entries are linked to the experiments table.

Appendix B

Figure B1 Weighted Spearman rank correlation coefficients between K_h at different tensions and climatic features, soil properties, land use and management factors and methodological details. In contrast to Fig. 10, only the 193 data entries from arable fields using a dry-to-wet sequence and the steady-state piecewise method (Reynolds and Elrick, 1991; Ankeny et al., 1991) were considered. Correlation coefficients are shown up to p values of 0.001.