Introduction

Soils generally consist of several and differently structured horizons and present direction-dependent properties due to the textural composition and this soil aggregate formation (Horn and Kutilek 2009; Alaoui et al. 2011). Its physical properties may vary from one location to another (heterogeneity), as well as from one direction to another at the same location (anisotropy). Even in the same layer, the heterogeneous arrangements of soil particles and aggregates due to anthropogenic and pedogenetic effects result in the discrepancy of pore continuity and tortuosity between the different directions (Dörner and Horn 2009; Beck-Broichsitter et al. 2020a). The directional behaviour (isotropy or anisotropy) governs not only hydrologic processes (e.g., infiltration, recharge) but also soil aeration and root growth (Beck-Broichsitter et al. 2020c; Zhang and Peng 2021).

The saturated hydraulic conductivity, Ks, and the corresponding soil water diffusivity are of primary importance in analysing rates and directions of water flow and solute transport in the vadose zone in agricultural and forest hillslope catchments (Weiler and McDonnell 2007; Dohnal et al. 2012; Horn et al. 2014; Dusek and Vogel 2018). It is important to note that analysis of soil hydraulic properties on the catchment scale are regularly expensive and time consuming. Thus, experiments on the small-scale (e.g., soil pits) for determining Ks values are preferably used for model-based upscaling (Weiler et al. 2003; Dusek and Vogel 2016). The method of upscaling is associated with uncertainties and limited significance in heterogenous catchments and hillslopes, although large-scale experiments represent an efficient way to obtain directional soil hydraulic properties functionally averaging local heterogeneities (e.g., Brooks et al. 2004; Pirastru et al. 2017). Other studies found that the occurrence of preferential flow in hillslope catchments are governed by small-scale soil structures (Carminati et al. 2007; Wiekenkamp et al. 2016). These local differences can be caused by tree roots, biopores, (Dohnal et al. 2012), shrinkage cracks and differently sized stones (Dusek and Vogel 2019), and tree- and machinery induced layering through compaction (Godefroid and Koedam 2010; Beck-Broichsitter et al. 2020b,c).

The X-ray computed tomography (μCT) is used to obtain 3D images and visualize the distribution and morphology of soil macropore networks and its characteristics (i.e., connectivity) for analysing local flow processes in macropores (i.e., Chandrasekhar et al. 2019; Leue et al. 2020; Schlüter et al. 2020). By contrast, additional methods as for example the scanning electron microscope (SEM) technique are more suitable for analysing the microstructure of soils in high-resolution including anisotropy and spatial variability of the permeability in unsaturated soils (Xu et al. 2020, 2021). Comparing both techniques, the μCT scanner was the method of choice with focus on the macropore network of undisturbed soil samples with a volume of 280 cm3.

The current study shows by means of a concrete example small-scale direction-dependent differences of a BwC horizon of a BwC horizon of a Dystric Cambisol under spruce forest formed on the weathered and fractured granite bedrock in the mountainous hillslopes Uhlirska catchment (Czech Republic). The idea of the study is to combine the soil capacity-related volume fraction of the pore classes, h(θ), with the intensity-related hydraulic conductivity function, K(θ), through the soil water diffusivity, D(θ), and their direction-dependent behaviour by an extended anisotropy factor, i.e., D(θ)*A. The aim was to identify the effect of small-scale soil structures in terms of cracks and differently sized stones affecting the soil water diffusivity between the macropore and fine pore size range for better understanding of the process affecting subsurface lateral flow in a forest hillslope. As a working hypothesis, a pronounced impact of macropores and wide coarse pore size range on D(θ) function was expected.

Materials and methods

Study site

The experimental headwater catchment Uhlirska is in the northern part of the Czech Republic at the Jizera Mountains (15′18″49″E, 50′49″26″N). The area of the catchment is around 1.78 km2 (Fig. 1). The mean annual temperature is 4.7 °C and the mean annual precipitation is 1380 mm. The elevation ranges from 779 m above the mean sea level at the closing profile to 886 m a.s.l. at the catchment divide; the average altitude is 822 m a.s.l. The catchment valley is formed by the Cerna Nisa stream and gentle slopes. The average length of the hillslopes is 450 m and the slope angle vary between 5 and 20% (Sanda and Cislerová 2009). The soils at hillslopes are shallow sandy loams classified as Cryptopodzols, Cambisols (Fig. 2), and Podzols (Sanda et al. 2014; Beck-Broichsitter et al. 2022). The catchment was deforested due to acid rain in the 1980s and 1990s and then reforested in 1995 by spruce monoculture (Picea abies). The catchment has been equipped with an extensive monitoring and sampling network (Sanda et al. 2014).

Fig. 1
figure 1

Scheme of Uhlířská catchment, Jizera Mountains (15′08″46.2″E, 50′49″58.5″N), Czech Republic. Modified after Sanda et al. (2014). Map of Europe (open source: Eurostat)

Full size image
Fig. 2
figure 2

Dystric Cambisol with AE/Bs/Bsw/BwC/C horizons located in the Jizera Mountains in northern Bohemia in the Czech Republic (Photo by M. Leue). Undisturbed soil cores in vertical (0°, y = x-axis) and horizontal (90°, z-axis) sampling direction

Full size image

Soil classification and basic soil characteristics

The hillslope soil is classified as Dystric Cambisol, formed on the weathered and fractured granite bedrock with AE/Bs/Bsw/BwC/C horizons (IUSS Working Group WRB, 2015). The disturbed and undisturbed soil samples were collected in September 2020 from one soil pit at the hillslope with a mean gradient of 14% (8°). Soil material of each horizon was also used to investigate the organic carbon (OC) content at 1200 °C and carbonate through Scheibler method. Soil texture (stones > 2 mm, sand: 0.063 mm to 2 mm; silt: 0.002 mm to 0.063 mm; clay: < 0.002 mm) was analysed by the combined sieve and pipette method with three repetitions each (Hartge and Horn 2016) and further classified by FAO (2006).

Soil physical properties

In total, 24 stainless steel rings (280 cm3, height: 5.5 cm; diameter: 8 cm) of undisturbed soil samples were used, 12 for vertical (0°, y = x-axis) and horizontal (90°, z-axis) sampling direction, respectively (Fig. 2). The samples were taken from the BwC horizon in 0.35–0.50 cm depth. At laboratory, the sample rings were first carefully wetted for determining the saturated hydraulic conductivity, Ks (m s−1), and in a second step, the water retention curve and the unsaturated hydraulic conductivity, K(h), were simultaneously determined. In a last step, the sample rings were dried for 24 h at 105 °C to determine the residual amount of water content, θr (cm3 cm−3), and the dry bulk density, ρ (g cm−3), respectively, while initial water content, θ was determined by total water loss related to core volume.

The KSAT device (Meter Group AG, Munich, Germany) was used to derive the Ks values based on the German method standards (DIN 18130-1 1998). The Ks values were determined in the falling head mode and fitted by an exponential function using the KSAT VIEW software version 1.5.0 (Meter Group AG, Munich, Germany):

$$h\left(t\right)={h}_{0}\mathrm{exp}\left(-K\mathrm{s}\bullet \frac{{A}_{\mathrm{bur}}}{{A}_{\mathrm{sample}}}\bullet \frac{1}{L}\bullet t\right)=a\mathrm{exp}\left(-b\bullet t\right)+c,$$
(1)

where h(t) is the pressure head (cm) at a certain time, t(s), h0 is the pressure head at t = 0, Abur is the cross-sectional area of the burette (4.524 cm2), Asample is the cross-sectional area of the soil sample (50 cm2), L is the length of the soil sample (5 cm), Ks is the saturated hydraulic conductivity (cm s−1), and a (cm), b (s−1), and c (cm) are the fitting parameters of the fitted exponential pressure head function.

The HYPROP measurement device (Meter Group AG, Munich, Germany) was used for simultaneous determination of the soil water retention curve and the hydraulic conductivity function (Schindler et al. 2010) using the extended evaporation method (Schindler 1980). The wetted sample rings were placed on the HYPROP base containing two vertically installed tensiometer at 1.25 cm and 3.75 cm distance from the bottom, both sealed at the bottom and placed on a balance (UMS GmbH Munich 2015). The sample mass changes (g) and pressure (hPa) in intervals of 10 min were automatically recorded by HYPROP-VIEW software version 1.4.0 (Meter Group AG, Munich, Germany). The samples dried out to about –600 hPa (lower tensiometer) and –850 hPa (upper tensiometer).

The h values at the upper, hu(t), and the lower, hl(t), tensiometer were used to determine the hydraulic gradient, im (cm cm−1), per time interval, ∆t = t0 – t1, along the horizontal distance, ∆z, between both tensiometer (Schindler 1980; Wendroth et al. 1993):

$${i}_{m}=\frac{1}{2}\left(\frac{{h}_{\mathrm{u}}\left({t}_{0}\right)-{h}_{l}\left({t}_{0}\right)}{\Delta z}+\frac{{h}_{u}\left({t}_{1}\right)-{h}_{\mathrm{l}}\left({t}_{1}\right)}{\Delta z}\right)-1.$$
(2)

The K(h) values were calculated according to the Darcy–Buckingham law (Schindler 1980; Wendroth et al. 1993):

$$K\left({h}^{*}\right)=\frac{\Delta m}{{\alpha A}_{s\mathrm{ample}}\cdot {\rho }_{\mathrm{H}2\mathrm{O}}\bullet \Delta t\cdot {\mathrm{i}}_{m}},$$
(3)

where h* (hPa) is the mean value of the pressure heads of the two tensiometer (i.e., h* = 0.5(hu + hl)) with ∆t = 10 min, ∆m is the sample mass difference (g) in the time interval, α is a flux factor (-), ρH2O is the density of water (~ 1 g cm3).

The data points of the soil water retention curve were calculated based on the water loss per volume of the sample and the corresponding geometric mean tension of the sample in the same time interval (UMS GmbH Munich 2015). Pressure chamber was used to determine moisture content at h value of -15,000 hPa.

The HYPROP-Fit software was used for calculation, evaluation, and fitting of the observed h(θ) and K(h) functions of the vertical and horizontal samples. The observed water retention data were fitted with the constraint unimodal (van Genuchten 1980) and bimodal van Genuchten model (Durner 1994):

$$\mathrm{Se}\left(h\right)=\frac{\theta -{\theta }_{r}}{{\theta }_{s}-{\theta }_{r}}=\left\{\begin{array}{c}{\sum }_{\mathrm{i}=1}^{2}{\mathrm{w}}_{\mathrm{i}}{\left[1+{\left({\alpha }_{i}h\right)}^{{n}_{i}}\right]}^{{-m}_{\mathrm{i}}} for \left|\mathrm{h}\right|<0\\ 1 for \left|\mathrm{h}\right|\ge 0\end{array}\right\},$$
(4)

where Se(h) is the effective saturation (-), θs is saturated and θr is residual water content, respectively, α (cm−1), n (-) and m (-), with m = 1—1/n, are empirical van Genuchten fitting parameters for both pore domains (index i), and wi is a pore-domain weighting factor (bimodal: w1 = 1—w2, unimodal: w2 = 0).

The unsaturated hydraulic conductivity, K(Se) was determined as function of the effective saturation with the combined Mualem–van Genuchten type analytical function and the bimodal retention function (Priesack and Durner 2006):

$$K\left(\mathrm{Se}\right)=K\mathrm{s}{\left(\sum_{\mathrm{i}=1}^{\mathrm{k}}{\mathrm{w}}_{\mathrm{i}}{\mathrm{Se}}_{\mathrm{i}}\right)}^{\uptau }{\left(\frac{\sum_{\mathrm{i}=1}^{\mathrm{k}}{\mathrm{w}}_{\mathrm{i}}{\mathrm{\alpha }}_{\mathrm{i}}\left[1-{\left({1-{\mathrm{Se}}_{\mathrm{i}}}^{1/{\mathrm{m}}_{\mathrm{i}}}\right)}^{{\mathrm{m}}_{\mathrm{i}}}\right]}{\sum_{\mathrm{i}=1}^{\mathrm{k}}{\mathrm{w}}_{\mathrm{i}}{\mathrm{\alpha }}_{\mathrm{i}}}\right),}^{2}$$
(5)

where integer k denotes the modality of the model (i.e., k = 2 for bimodal), wi is a weighting factor for the sub-curves of each pore domain (0 < wi < 1 and ∑wi = 1), and τ (-) is the tortuosity.

In a first step, the water retention function was fitted as well as the derived fitting parameters αi, ni, mi, which were then used in a second step for fitting the K(Se) functions.

The most appropriate soil water retention model is selected by the second-order Akaike Information Criterion (Hurvich and Tsai 1989) was used. The error between fitted and observed K(Se) values is automatically determined by the root mean square error (RMSE) through HYPROP-Fit software.

The pore size distribution (PSD) was calculated by the slope of the water retention curve ΔSe/Δh. The h value of − 10 hPa (about 9.81 cm) was assumed as boundary between the macropores MaP (h1 ≥ − 10 hPa), structural pores: wide coarse pores wCP (− 10 > h2 ≥ − 60 hPa), and textural pores: narrow coarse pores nCP (− 60 > h3 ≥ -300 hPa), medium pores MP (− 300 > h4 ≥ -15,000 hPa), and fine pores FP (h5 < − 15,000 hPa) assumed from the peaks in the PSD at an equivalent pore diameter (Beck-Broichsitter et al. 2020a,b). The pore diameter for the different pore size classes were calculated as (Hartge and Horn 2016)

$${r}_{\mathrm{i}}=\frac{0.149}{\left|{\mathrm{h}}_{\mathrm{i}}\right|},$$
(6)

where ri is the pore radii (cm) and h is the pressure head (cm), assuming 1 cm is proportional to − 1 hPa, for the ith pore size class: MaP (h1 ≥ -10 hPa), wCP (− 10 > h2 ≥ − 60 hPa), nCP (− 60 > h3 ≥ − 300 hPa), MP (− 300 > h4 ≥ -15,000 hPa), and FP (h5 < -15,000 hPa). Using the Laplace equation for a capillary rise, these five pore classes can be defined with pore diameter d as: MaP with d > 0.3 mm, wCP with 0.3 ≥ d > 0.05 mm, nCP with 0.05 ≥ d > 0.01 mm, MP with 0.01 ≥ d > 0.0002 mm, and FP with d ≤ 0.0002 mm.

The anisotropy ratio (AR) is defined as ratio of K(Se)-values obtained in horizontal (hor) and vertical (ver) direction (Beck-Broichsitter et al. 2020a, b):

$$\mathrm{AR}=\frac{{K(\mathrm{Se})}_{\mathrm{hor}}}{{K(\mathrm{Se})}_{\mathrm{ver}}},$$
(7)

where the subscripts ver and hor indicate K(Se) values obtained from cores sampled in vertical and horizontal directions, respectively.

The soil water diffusivity, D(θ), is used for linking the intensity parameter with the capacity-based volume fraction of ith pore size classes with a pore size-related fraction of the capillary potential (Beck-Broichsitter et al. 2020a). The D(θ) (cm2 d−1) corresponds to the product of the K(θ) and the differential inverse water capacity, dh/dθ (Hillel 1998):

$$D\left(\theta \right)=K\left(\theta \right)\frac{\mathrm{dh}}{\mathrm{d\theta }},$$
(8)

The effective D(θ) values for the ith pore size classes were calculated as (Beck-Broichsitter et al. 2020a)

$${D}_{\mathrm{i}}\left({\theta }_{\mathrm{i}}\right)={\int }_{\uptheta +}^{\uptheta -}K\left({\uptheta }_{\mathrm{i}}\right)\frac{\mathrm{dh}}{\mathrm{d}{\uptheta }_{\mathrm{i}}}\mathrm{d\theta },$$
(9)

where θ+(h+) is the upper and θ─(h─) is the lower value of each of the five h-ranges defined above.

A diffusivity-dependent effective value of the anisotropy ratio ARD(θ) (cm2 d−1), per pore size class i was calculated as (Beck-Broichsitter et al. 2020a)

$${{\mathrm{AR}}_{\mathrm{D}(\uptheta )}=\mathrm{D}}_{\mathrm{i}}({\uptheta }_{\mathrm{i}})\underset{{\uptheta }^{+}}{\overset{{\uptheta }^{-}}{\int }}\frac{1}{{\mathrm{AR}}_{\mathrm{i}}}\mathrm{d\theta }.$$
(10)

The diffusivity-dependent anisotropy is a direct link between intensity and capacity parameters and improves the information about the impact of the pore size classes on Ks and finally D(θ) values (Beck-Broichsitter et al. 2020a).

X‐ray computer tomography

Internal structure of the two field moist undisturbed soil samples in stainless steel cores (280 cm3, height: 5.5 cm; diameter: 8 cm) (withdrawn in vertical and horizontal direction) with volume size of 1741*1854*1621 voxel was visualized by X-ray computer tomography (CT). The CT images were obtained using a GE v|tome|x M 300 μCT scanner at the Technical University in Brno (Czech Republic). The X-ray photon energy was set to 220 keV with a current of 300 μA (GE Sensing & Inspection Technologies GmbH, Wunstorf, Germany). A 1-mm Cu pre-filter was used to minimize beam hardening. Each soil sample was scanned in two height steps resulting in two overlapping 3D image data sets that were subsequently reconstructed to obtain a single 3D core volume. For reconstruction, the software datos|x, version 1.5.0.15 (GE Sensing Technologies) was used, which is based on a modified Feldkamp algorithm (Feldkamp et al. 1984). A spatial resolution in terms of an isotropic voxel edge length of 55 μm was achieved for both samples.

The image pre- and post-processing was done with the SoilJ software (Koestel, 2018) implemented as plugin for the open-source software ImageJ 1.x (Schneider et al. 2012). About 2 mm of the outermost layer of the soil samples were clipped out by image grayscaling to remove the cylinder wall to calculate the bulk density of the region of interest (ROI). This information can also be used for SoilJ’s beam-hardening correction routine to remove scattering artifacts from images of steel rings (Hansson et al. 2017).

SoilJ offers image segmentation options as well as the possibility of analysing joint histograms of several calibrated images (Koestel 2018) and for the binarization of the gray value images, a global threshold was calculated for each of the 403 slices per soil core (ROIvertical) and 407 slices per soil core (ROIhorizontal), using Otsu’s method (Otsu 1979). Histograms of objects volumes (e.g., pores, soil matrix, stones) were used to determine thresholds between the pore volume classes (Pagenkemper et al. 2013; Leue et al. 2020); objects with diameter of around 1.134 pixel were denoted as macropores (d > 0.3 mm).

The global morphological measures of the selected ROI (critical diameter and thickness of pores, macropore volume) were calculated in addition to the Minkowski functionals. For this purpose, SoilJ makes use of the algorithms included in the MorphoLibJ (Legland et al. 2016) and BoneJ (Doube et al. 2010) packages.

To quantify the anisotropy and, therefore, the directionality of heterogenous structured materials, tensorial structural characteristics are needed to determine the degree of anisotropy and the preferred orientation. The mean intercept length (MIL) tensor is the most used technique for it (Odgaard 1997; Klatt et al. 2017) and the degree of the CT-based anisotropy (ACT) is defined as follows (Doube et al. 2010):

$$A_{{{\text{CT}}}} = 1 - \frac{{\text{Long axis eigenvalue}}}{{\text{Short axis eigenvalue}}},$$
(13)

where ACT = 0 means isotropic and ACT = 1 anisotropic condition, and the longest axis has the smallest eigenvalue.

The Paraview software version 5.9.2-RC2 (e.g., Ahrens et al. 2005) was used for visualisation of the macropore volumes (MaPCT: d > 1.134 pixel) of the vertical and horizontal sampled soil cores.

Results and discussion

Soil physical and chemical properties

The soil texture is classified as loam for the AE/Bs/Bsw horizons and sandy loam for the BwC horizon with clay contents between 110 and 120 g kg−1 and stone contents between 87 and 240 g kg−1 (Table 1). The pH values vary between 3.9 and 4.5, and the organic carbon content ranges from 35 in the BwC horizon to 65 g kg−1 in the Bs horizon. The ρ values increase from 0.65 to 1.25 g cm−3 with increasing soil depth.

Table 1 Soil physical and chemical properties of the AE/Bs/Bsw/BwC/C horizons of the Dystric Cambisol; organic carbon (OC) content, pHCaCl2, and soil texture, dry bulk density (ρ); mean values and standard deviations (symbol ±) of three repetitions each, five for ρ
Full size table

The results in Table 2 show higher Ks values with 554 cm d−1 for the horizontal direction compared to the vertical one with 123 cm d−1 with ρ values of 1.23 g cm−3, respectively.

Table 2 Total porosity, TP (cm3 cm−3) and volume fractions of pore classes for macropores MaP (h1 ≥ − 10 hPa), b) wide coarse pores wCP (− 10 > h2 ≥ − 60 hPa), c) narrow coarse pores nCP (− 60 > h3 ≥ − 300 hPa), d) medium pores MP (− 300 > h4 ≥ -15,000 hPa), and e) fine pores FP (h5 < -15,000 hPa)
Full size table

The volume fractions of pore size tend to be higher for the horizontal than vertical collected samples, especially the volume fraction of macropores. MaP is more pronounced for the horizontal direction with 0.079 cm3 cm−3 than the vertical one with 0.038 cm3 cm−3. With respect to the total porosity, the macropore volume is around 7.4% in vertical and 13.1% in horizontal direction (Table 2).

The RMSE values in Table 3 vary between 0.0012 and 0.002 cm3 cm−3 for θ and between 0.1188 and 1.148 cm d−1 for Log K(Se). The AICc values for the bimodal model are smaller (larger in absolute value) in both directions than for the unimodal model. Thus, the bimodal model is found to be provide more appropriate (Table 3).

Table 3 Root mean square error (RMSE) for mean deviation between the fitted and observed water content θ, and the unsaturated hydraulic conductivity Log K(Se)
Full size table

The results in Table 4 indicate that the fraction of the second pore domain, w2, representing the smaller pore size is notable higher in vertical direction with 0.724 than in horizontal direction with a more dominating larger pore size (1–w2) with 0.809.

Table 4 Hydraulic parameter of the soil water retention data with bimodal van Genuchten model (Durner 1994) and bimodal unsaturated hydraulic conductivity K(Se) (Priesack and Durner 2006) in vertical (ver) and horizontal (hor) direction; mean and standard deviations (± symbol) from 12 replicates
Full size table

The slopes of the soil water retention curves, ΔSe/Δh, per interval in vertical and horizontal direction shown in Fig. 3c, indicate a structural peak at − 10 hPa. A more pronounced second matric peak at about − 1000 hPa is seen for the vertical direction (Fig. 3c), while the second matric peak for the horizontal direction at about − 300 hPa is less pronounced (Fig. 3d). Overlooking the bi-peak structure, the higher the amount of macropores and medium pores, the more pronounced is the structural peak as also mentioned in Ding et al. (2016) and Beck-Broichsitter et al. (2020b).

Fig. 3
figure 3

Soil water retention curves approximated parameter fits from the mean of 5 data points with HYPROP-Fit software of the average and two representative soil water retention data sets with bimodal (w1 = 1–w2) constrained (m = 1–n-1) van Genuchten model (Durner 1994); corresponding relative pore volume fractions (1c,d) at the measured pressure head intervals in terms of the mean slope of the effective saturation (ΔSe/Δh) per interval in vertical and horizontal direction

Full size image

The results indicate that the used HYPROP measurement system is an appropriate tool for soils with bimodal pore size distribution (Durner 1994) considering two different flow domains: the macropore system and the soil matrix (Weninger et al. 2018). The direction-dependent differences in the pore size distribution as shown in Fig. 3c and d differ from the theory of isotropic capacity parameter as also disproved by Horn et al. (2014). However, differences in the volume fraction of the pore classes and, therefore, the peak structure can be related to the small-scale heterogeneity in soil profiles (Horn et al. 2019; Beck-Broichsitter et al. 2020a, b). On the other hand, capturing a larger volume fraction of the pore classes is greater for vertically as compared to horizontal-oriented soil sampling (Beck-Broichsitter et al. 2020a).

Unsaturated hydraulic conductivity and anisotropy ratio

The functions of K(Se), fitted with the bimodal model are shown in Fig. 4. The results show weak direction-dependent differences of K(Se) between the two directions of sampling. The values of K(Se) are slightly larger in horizontal than in vertical direction for saturated conditions (Se ranging from 0.8 to 1).

Fig. 4
figure 4

Functions of K(Se) of the average and two representative K(Se) functions with bimodal retention function in vertical and horizontal direction

Full size image

The results in Table 2 indicate higher mean Ks values in horizontal (554 cm d−1) than in vertical direction (123 cm d−1) and in Fig. 5, the anisotropy ratio AR is shown as a function of pressure head. In the range of the macropores (h ≥ − 10 hPa), the AR values decrease from 4.6 to 0.3. With the emptying of the soil pores, the AR values increase and peak at about − 60 hPa. The values of AR increase up to 3.5 in the unsaturated conditions.

Fig. 5
figure 5

Average anisotropy ratio (K(Se)hor/K(Se)ver) plus maximum ARmax (highest K(Se)hor/K(Se)ver) and minimum ARmin (lowest K(Se)hor/K(Se)ver) from 12 replicates each as function of the pressure head

Full size image

The anisotropy ratio of 4.5 at saturated conditions suggests an increased potential for lateral flow above a less permeable bedrock. Predominantly horizontal-oriented pore network was created by plant roots and root fragments (Sanda and Cislerová 2009; Dohnal et al. 2012), and differently sized stones and rock fragments (Dusek and Vogel 2019). This assumption is further underlined by a macropore volume of 8.2% in vertical and 14.8% in horizontal direction (Table 2).

As recently mentioned, the anisotropy of hydraulic conductivity is often neglected, because it can only partially explain the complex variably saturated response dynamics in layered hillslope soils (Mirus 2015).

Soil water diffusivity

The functions of soil water diffusivity D(θ) for the volume fraction of the pore classes show that the function of D(θ) has a positive slope for the macropores and wide coarse pores with an r2 between 0.66 and 0.84 for both sampling directions (Fig. 6). On the other hand, the narrow coarse pores, medium pores, and fine pores show less pronounced positive and also negative linear regressions with r2 between 0.29 and 0.44.

Fig. 6
figure 6

Soil water diffusivity, D(θ), as a function of volume fraction of five pore classes: (a) macropores MaP (h1 ≥ − 10 hPa), (b) wide coarse pores wCP (− 10 > h2 ≥ − 60 hPa), (c) narrow coarse pores nCP (− 60 > h3 ≥ -300 hPa), (d) medium pores MP (− 300 > h4 ≥ -15,000 hPa), and (e) fine pores FP (h5 < -15,000 hPa). The dashed lines indicate the confidence intervals of 95% of the relations between the D(θ)-values and the pore size distribution (pore size) fitted to 12 data points in vertical (ver) and horizontal (hor) direction, respectively

Full size image

In Fig. 7, the D(θ)-weighted anisotropy ratios ARD(θ) defined in Eq. 10 are depicted. The results show a positive linear relationship for the MaP (r2 of 0.89–0.92) and wCP (r2 of 0.49–0.61), while nCP, MP, and FP indicate a negative linear relationship (r2 of 0.28–0.52). Thus, the anisotropy ratios ARD(θ) strongly increases with the volume fraction of macropores and wide coarse pores, respectively.

Fig. 7
figure 7

D(θ)-weighted average anisotropy ratio ARD(θ) as a function of the volume fraction of five pore classes: (a) macropores MaP (h1 ≥ − 10 hPa), (b) wide coarse pores wCP (− 10 > h2 ≥ − 60 hPa), (c) narrow coarse pores nCP (− 60 > h3 ≥ − 300 hPa), (d) medium pores MP (− 300 > h4 ≥ -15,000 hPa), and (e) fine pores FP (h5 < -15,000 hPa). The dashed lines indicate the confidence intervals of 95% of the relations between ARD(θ) and the volume fraction of pore classes fitted to 12 data points in vertical (ver) and horizontal (hor) direction, respectively

Full size image

The linear regression between D(θ) and the diffusivity-dependent effective value of the anisotropy ratio, ARD(θ), for the narrow coarse to fine pore range shows r2 < 0.52 (Fig. 7). Thus, intensity parameters are also functions of pore continuity and connectivity and not directly of capacities (Horn et al. 2009, 2014) as exemplary shown in Fig. 8.

Fig. 8
figure 8

Vertical cross section of automatically detected soil column outlines (red) of CT-based images for the vertical sampled soil core (a) and the horizontal sampled soil core (b) for further validation purposes

Full size image

The K(Se)- and, therefore, the D(θ) values are mostly related to the Ks values (Germer and Braun 2015) and the volume fraction of the macropores and wide coarse pores (Beck-Broichsitter et al. 2020a, b). The higher van Genuchten parameter n1 and n2 in vertical direction (1.766 and 1.504) than n1 and n2 in horizontal direction (1.402 and 1.225, see Table 4) indicate a lower water capillary adsorption capacity for the macropores and the coarse to medium pores in vertical compared to horizontal direction. As a result, a more pronounced drop in conducting water when drained is observed for vertical than horizontal direction (Alaoui et al. 2011; Beck-Broichsitter et al. 2020a).

It is suggested that the BwC horizon represents a critical part of the hillslope soil profile with hydraulic characteristics favouring generation of subsurface lateral preferential flow at the mountainous Uhlirska catchment, especially during snow melt events and heavy rain storms that provide nearly saturated soil conditions at the base of the soil profile (Sanda and Cislerova 2009; Sanda et al. 2014). In addition, the hillslope provides the slope gradient that is needed for lateral subsurface water flow (Widomski et al. 2015; Beck-Broichsitter et al. 2018). This assumption is underlined by a study of Wiekenkamp et al. (2016), which confirmed that the spatial and temporal occurrence of subsurface lateral flow is triggered by the amount and intensity of precipitation (Weiler and McDonnell 2004).

It is mentioned by Brooks et al. (2004) that small-scale Ks values do not represent large scale Ks value due to a general lack in a necessary number of samples and also not reproduce larger pore networks. By contrast, machine learning algorithms can be used to develop high-performance pedotransfer functions regarding the sensitivity of Ks values to the soil structure (Araya and Ghezzehei 2019). Thus, large-scale Ks values from in situ experiments can be used for training the algorithms to improve the predicting function of small-scale Ks values and, therefore, the model-based upscaling performance.

The D(θ) as direct link between K(Se) values and the pore size classes improve the validity for describing pore networks on small scales (Beck-Broichsitter et al. 2020a) compared to standalone Ks values that can be hardly predicted by pedotransfer functions (Mirus 2015; Pirastru et al. 2017).

CT image analysis

The process of image segmentation and the pore space analysis using SoilJ (Koestel et al. 2018) show a CT-based macropore volume (MaPCT) of 0.028 cm3 cm−3 for the vertical and 0.047 cm3 cm−3 for the horizontal oriented soil cores (Table 5). The MaPCT values are slightly lower than the mean MaP values in Table 2 with 0.038 cm3 cm−3 and 0.079 cm3 cm−3, respectively. The discrepancy was probably caused by different measurement approaches (e.g., Pagenkemper et al. 2013) and a threshold-based underestimation of macropore volume (Koestel 2018; Leue et al. 2020). The CT technique and the used SoilJ software (Koestel et al. 2018) enabled a rapid extraction and quantification of structural features. The critical pore diameter is the bottleneck in the connection from top to bottom and varies between 9.28 vx (vertical) and 9.15 vx (horizontal).

Table 5 Results of the CT image analysis for the vertical (ver) and horizontal (hor) sampled soil core using SoilJ (Koestel 2018)
Full size table

The ACT values considering the MaPCT are larger for the horizontal sampled soil core than for the vertical, and both are classified as predominantly isotropic, while the predominant orientation of macropores is x-axis (= y-axis) for the vertical and z-axis for the horizontal oriented soil core and, therefore, parallel to the sampling direction (90°). (Fig. 9c, d).

Fig. 9
figure 9

2D images (8bit) of the top view (a, b) for vertical (0°) and horizontal (90°) sampled soil cores (height: 5.5 cm; diameter: 8 cm). 3D images of the vertical (c) and horizontal (d) samples. The 3D views of the vertical (c) and horizontal (d) oriented macropore volume (d > 1.134 pixel) with an image resolution of 1407*1407 pixel in x/y extent and 395 pixels in z extent, respectively. The 3D images were created with the visualisation software Paraview (e.g., Ahrens et al. 2005). Color’s RGB values > 200 indicate pores with diameter > 1 mm (~ 3.8 pixel)

Full size image

The results in Fig. 9 show that stones and rock fragments form an abandoned and continuous pore network from µm- to mm-scale. It is obviously pore structure which differs from earthworm- and root-derived macropores (e.g., Pagenkemper et al. 2013; Schlüter et al. 2020). Its connectivity depends on the layering and the size of stones and fragments, and even stones can have some porosity, based on the way they were weathered (e.g., Naseri et al. 2019). The predominant orientation of macropores in z-axis of the horizontal sampled soil core underlines the findings that weathered bedrock strata can initiate lateral preferential flow especially during intense rainfall events (Dusek and Vogel, 2014). However, further analysis is necessary using the derived results for upscaling of small-scale macropore network information upon the catchment scale as also proposed by Schlüter et al. (2020).

Anisotropic hydraulic conditions and slope stability

Considering the term soil structure again, anisotropic soil hydraulic properties may help in evaluating slope stability at the hillslope scale in relation to the rainfall characteristics and environmental conditions. Especially soil structure dependent permeability contrasts among soil horizons can impact infiltration processes and the pore water pressure distribution (Zhao and Zhang 2014) resulting in locally diverted flow dynamics (Formetta and Capparelli 2019; Fusco et al. 2021). Soil deformation due to mechanical loads can lead to a distinct horizontal anisotropy of the pore system in the (sub)soil because of anthropogenic induced plate structure formation. Accompanying positive pore water pressures can result in a decrease of the shear strength or an increase in driving forces causing slope instability (Yeh et al. 2015; Hartge and Horn 2016; Fusco et al. 2021). Thus, it becomes obvious, that shear parameters are tensors which depend on the aggregation besides all further internal properties (Dörner and Horn 2009).

However, additional research is needed to combine anisotropy of soil hydraulic and mechanical properties for evaluation slope stability under various soil and vegetation types.

Conclusion

The results of this study indicate direction-dependent differences in the soil hydraulic conductivity as function of effective saturation, K(Se). The soil water diffusivity, D(θ), as link between the intensity and capacity parameter was found to strongly depend on the pore size distribution. Thus, by weighting the anisotropy ratio with the soil water diffusivity, both the capacity and the intensity parameters were combined. The analysed X-ray CT images indicate that more pronounced cracks were observed in the soil samples taken in horizontal than in vertical direction in the BwC horizon of hillslope profile.

The anisotropy in the range of the soil macropores was determined, resulting in increased value of hydraulic conductivity in the horizontal direction. Similar trend was found for the medium pores. In the current study, the anisotropy ratios above 1 in the volume fraction of the macropores close to saturation indicate exemplary the tendency of the BwC horizon of the Dystric Cambisol for supporting the subsurface lateral water flow of the mountainous hillslopes Uhlirska catchment. Based on the study results, the refined information on soil hydraulic characteristics have the potential to improve runoff predictions at the hillslope scale. The X-ray computer tomography-based anisotropy ratio is larger for the horizontal sampled soil core than for the vertical, and both are classified as predominantly isotropic. However, the predominant orientation of macropores in z-axis of the horizontally sampled soil core and, therefore, parallel to the sampling direction (90°) underlines the findings that weathered bedrock strata can initiate lateral preferential flow.

In future, the small-scale direction-dependent differences in the soil hydraulic capacity and intensity parameter will be used for model-based upscaling for better understanding of preferential flow at the catchment scale.