Routing stemflow water through the soil via preferential flow: a dual-labelling approach with artificial tracers

Stemflow and its belowground funnelling along roots and macropores may play an important role in the soil moisture redistribution in forest environments. In this study, a stemflow experiment on Pinus sylvestris L. (Scots pine) used artificial tracers to view and quantify preferential flow after stemflow infiltration into the soil. A total of 41 L of water labelled with enriched deuterium and brilliant blue FCF were applied at a flow rate of 7 L h⁻¹ to the stem of a pine tree, which corresponds to the stemflow caused by about 50 mm of rainfall. Time domain reflectometry (TDR) probes were installed around the tree trunk to measure the high-resolution volumetric water content. A total of 1 d after the stemflow discharge, soil pits were dug in the different cardinal directions and at varying distances from the tree. Photographs were taken for imaging analysis to quantify preferential flow metrics. Soil samples were taken from the different profiles to analyse the dye concentrations and isotopic compositions. We found that stemflow infiltrated through an annulus-shaped area around the tree base. We observed a heterogenous spatiotemporal soil moisture response to stemflow and the occurrence of shallow perched water tables around the tree trunk. Dye staining demonstrated that stemflow infiltrated primarily along the surface of coarse roots and through macropores. The dye coverage was less extensive close to the soil surface and increased with depth and with proximity to the tree trunk. Lateral flow was also observed, mainly in the shallow soil layers. Our analyses demonstrate the prevalence of preferential flow. Deuterium and brilliant blue FCF concentrations were significantly correlated. The tracer concentrations decreased with increasing distance from the tree trunk, indicating dilution and mixing with residual soil water. Macropores, coarse roots (living or decayed) and perched water tables produced a complex network regulating the preferential flow. Our results suggest that stemflow affects soil moisture distribution, and thus likely also groundwater recharge and surface runoff. Our study provides insights into the soil hydrological processes that are regulated by stemflow belowground funnelling and improves our understanding of forest–water interactions.

1 Introduction

Rainfall interception by forest canopies prevents rainfall from immediately reaching the soil, redistributing it in the form of throughfall and stemflow. Throughfall refers to the rainwater that may or may not have contact with vegetative surfaces and falls or drips to the forest ground, whereas stemflow is the proportion of rainwater that is channelled towards the bole and eventually flows downward and reaches the forest floor around the stem base. Stemflow contributes higher amounts of water and solutes in near-stem soils than the throughfall contribution in the rest of the forest (Carlyle-Moses et al., 2018). Stemflow varies widely across vegetation types and climatic regions (Levia and Frost, 2003; Llorens and Domingo, 2007; Yue et al., 2021). Recent research has unveiled biotic and abiotic factors that regulate stemflow production before it is funnelled belowground (Levia and Germer, 2015; Cayuela et al., 2018; Zhang et al., 2021). Although stemflow is not a major proportion of the overall catchment water balance, it can saturate near-stem soil and may enhance overland or subsurface matrix flow; contribute to soil water replenishment and groundwater recharge; produce soil erosion (Levia and Germer, 2015), transport chemicals and colloids (Van Stan et al., 2021); and influence root water uptake (Hildebrandt, 2020). Due to its importance in soil hydrology, biogeochemistry and ecohydrology, stemflow has become an active topic of research in recent years (Levia et al., 2011; Tucker et al., 2020; Van Stan et al., 2020; Magyar et al., 2021; Pinos et al., 2021; Tonello et al., 2021).

Root-induced bypass flow, introduced by Johnson and Lehmann (2006) as the “double-funnelling of trees”, refers to the process whereby an aboveground funnelling occurs when rainfall concentrates into stemflow, followed by a belowground funnelling process in which stemflow infiltrates into the soil along tree roots and macropores. Once stemflow occurs, its underground infiltration is governed by a complex interplay of soil physical properties, terrain slope and root architecture. While the concept of the double-funnelling phenomenon is well known, clarification of the underlying mechanisms is more challenging. Several experimental techniques, both invasive and non-invasive, have been used to study the stemflow belowground funnelling process of trees. The most common approaches are the monitoring of water fluxes by soil moisture sensors, the viewing of preferential flow paths by dye tracers or a combination of the two (Liang et al., 2011; Schwärzel et al., 2012; Spencer and van Meerveld, 2016; Gonzalez-Ollauri et al., 2020; Tischer et al., 2020). Recently, non-invasive geophysical techniques, such as ground-penetrating radar (GPR) and electrical resistivity tomography (ERT), have also been used to study belowground funnelling processes (Guo et al., 2020; Di Prima et al., 2022). The majority of these stemflow funnelling studies (of various deciduous and evergreen trees) found that stemflow infiltration leads to vertical and horizontal preferential flow along the surface of coarse roots, although matrix stemflow infiltration has been reported as well (Gonzalez-Ollauri et al., 2020).

Despite their potential for the study of water movement in the unsaturated zone, geophysical techniques often cannot resolve flow pathways at a small-scale resolution, which limits their usefulness in stemflow infiltration research (Fan et al., 2020). On the other hand, dye tracing for staining flow pathways has become an established way to show preferential flow in soils (Flury and Wai, 2003). The results obtained from staining experiments clearly illustrate the complicated patterns of water movement in the soil at a high spatial resolution (Ghodrati and Jury, 1990; Flury et al., 1994; Weiler and Flühler, 2004). However, the disadvantage of staining is that it is a destructive technique. Moreover, as dyes are adsorbed to some extent by the soil matrix (e.g. brilliant blue FCF; Ketelsen and Meyer-Windel, 1999; German-Heins and Flury, 2000), there is a certain degree of underestimation of the flow paths. Therefore, new approaches are needed to unravel the patterns of stemflow belowground infiltration. Stable isotopes, such as deuterium (²H), have been widely use in hydrology as conservative tracers (Kendall and McDonnell, 1998). In combination with dye tracers, they can be a powerful tool for unravelling stemflow belowground funnelling phenomena.

In this study, we developed an approach to characterize the stemflow belowground funnelling process by using high-resolution hydrometric monitoring and dual-tracer labelling (²H and dye). The main goals of our study were (1) to assess stemflow infiltration quantitatively by means of dye concentrations, isotopic compositions and hydrometric data (soil water content) and (2) to view the spatial distribution of preferential flow pathways within the soil by dye staining. The experiment was conducted in a Mediterranean ecosystem with pine forest cover, which had not been extensively investigated for stemflow belowground funnelling. This study is a contribution towards the mechanistic understanding of water infiltration patterns and water–root interactions in forested catchments.

2 Material and methods

2.1 Study site

The research was conducted at the Vallcebre research catchments in Catalonia, Spain, established as hydrological and ecohydrological research sites in 1988 (Llorens et al., 2018). The catchments are located 100 km north of Barcelona, in the southeastern part of the Pyrenees (42^∘12^′ N, 1^∘49^′ E). The climate is humid Mediterranean and characterized by a mean annual air temperature of 9.2 ^∘C and a mean annual precipitation of 856 mm (for the 1999–2018 period). Precipitation is seasonal, with spring and autumn being the wettest seasons and summer and winter the driest. Nevertheless, summer convective storms also provide significant precipitation. Reference evapotranspiration shows a seasonal pattern with a maximum in summer of up to 6.9 mm d⁻¹ and a mean annual rate of 823 mm. Soils developed over the mudstone lithology have a silty loam and silty clay loam texture and are characterized by a rapid decrease in their hydraulic conductivity at depth (Rubio et al., 2008).

Figure 1 Tracer experiment set-up. (a) Photographs of the artificial stemflow set-up and soil instrumentation before and after (last picture) the tracer experiment. (b) Top view of the experimental area showing microtopography and the location of instruments around the selected tree. The abbreviations used in the figure are as follows: swc – vertical 30 cm soil water content probe, swc_p – profile of four vertical soil water content probes, wt – mini-piezometer (30 cm deep), N – north, E – east, S – south and W – west. The number in the label of each location refers to the distance (cm) from the tree trunk.

Nowadays, Scots pine (Pinus sylvestris L.) forests, which arose through afforestation of old agricultural terraces, are the main land cover of the catchment, but small fragments of original oak (Quercus pubescens Willd.) forests also occur (Poyatos et al., 2003). An experimental Scots pine forest stand (named Cal Rotes) was delineated within the Can Vila catchment to investigate forest–water interactions. The stand density is 1189 trees per hectare and the basal area is 45.1 m² ha⁻¹. One mature pine tree (with a diameter at breast height of 27.3 cm and a basal area of 585.3 cm²), which had been previously instrumented for stemflow measurement during a 6-month period (Cayuela et al., 2018), was selected for our stemflow infiltration experiment. This tree is representative of the plot's diameter at breast height (DBH) distribution and is situated on a flat area on an old abandoned agricultural terrace. Further, it has sufficient space around it for the placement of all of the hydrometric equipment.

2.2 Artificial stemflow experimental set-up

To simulate stemflow (Fig. 1), the approach of Llorens et al. (2022) was adopted in this study. A flexible PVC plastic tube (internal diameter of 10 mm) was placed around the experimental tree at breast height (1.3 m). As the trunks of pine trees are very rough, the outer bark of the tree was removed to assist the assembly of the tube. The plastic tube was perforated on the inside with 2 mm holes every 4–5 cm, such that the holes were in contact with the trunk surface. To ensure uniform wetting of the stem and to prevent splashing of water, cloths were placed between the bole and the tube. The tube was connected by a Y-shaped tubing connector to a water tank. A stopcock was installed upstream of the Y-shaped tubing connector to regulate and control the flow rate.

Our system was designed to discharge water from a tank that was located 9 m above in an upper terrace (∼35 m away from the trunk), creating a gravity-fed water system. The volumetric soil water content (SWC, cm³ cm⁻³ or %) was monitored by sixteen 30 cm long time domain reflectometry (TDR) probes (CS615 probes, Campbell Scientific, Logan, UT, USA). The probes were placed vertically in two circular rings, 10 and 30 cm from the trunk every 45^∘, i.e. facing the cardinal and intercardinal directions (Fig. 1). When a physical obstacle was encountered (e.g. coarse roots or rocks) that did not allow the installation of the probe, the probe was placed as close as possible to the intended position. The TDR probes were measured every 20 s and averaged at 1 min intervals by a datalogger (Data Taker DT85, Thermo Fisher Scientific Inc., Waltham, MA, USA). The TDR probes had previously been calibrated with long-term manual TDR measurements performed with a Tektronix 1502-C cable tester (Tektronix Inc., Beaverton, OR, USA) at the study plot. In addition, in each main cardinal direction (N, S, E, and W), the SWC was monitored by four TDR probes of different lengths (CS605 probes, Campbell Scientific, Logan, UT, USA) that had been installed vertically at depths of 0–10, 0–20, 30–40 and 30–60 cm (Fig. 1) 2 years prior to the experiment. For their installation, holes with a diameter equal to the distance between the TDR probe external rods were made manually with a soil auger, and the TDR probe was stuck into the undisturbed soil once the hole had the selected depth. The probes were measured manually every 30 min with a Tektronix 1502-C cable tester.

Two piezometers were placed at 1.7 and 2.0 m depth at a distance 1.6 m from the tree bole in opposite directions (northwest and southeast, respectively) in order to monitor and sample groundwater. Water levels were measured every 20 min with dataloggers (Micro-Diver, Van Essen Instruments, Delft, South Holland, the Netherlands). Finally, twelve 30 cm deep mini-piezometers were placed at different distances (between 20 and 50 cm from the bole of the tree) in the cardinal and intercardinal directions to sample the soil water in the shallow layer (Fig. 1). When a physical obstacle was encountered, the mini-piezometer was placed as close as possible to the intended position. The water level was measured about every 30 min, and water samples were collected for all of the piezometers at 1 h intervals throughout the experiment.

Using the tree location as the centre point, a level mesh (4.8 m × 4.8 m with a grid size of 0.2 m) was placed above the ground. As depicted in Fig. 1, the mesh was used as a coordinate system to georeference the installed instruments and to characterize the local topography around the tree. The vertical distance (z) from the mesh intersection to the actual ground surface was measured, using the mesh level as the reference plane. All z values were corrected to represent the surrounding topography in relation to the tree base (z=0), and a digital elevation model (0.2 m resolution) was interpolated from these data, with the tree base as the reference elevation (Surfer v.18, Golden Software, Golden, CO, USA).

Unlike other stemflow infiltration studies (e.g. Schwärzel et al., 2012; Llorens et al., 2022), forest litter surrounding the tree bole was not removed prior to the experiment. The litter, which was 4.7±2.4 cm thick according to Molina et al. (2019), was left in place in order to better replicate natural stemflow infiltration conditions.

The stemflow experiment was conducted on 26 May 2021 when the soil had a moderate moisture condition (a mean prior volumetric soil water content of 32 %). During the previous 20 d, only 15 mm of rain was recorded, of which 6.4 mm occurred on 22 May. We applied 41 L of stemflow water labelled with enriched deuterium and brilliant blue FCF dye at a flow rate of 7 L h⁻¹. This amount and rate of stemflow correspond to that commonly generated by a 50 mm rainfall event for this particular tree (Cayuela et al., 2018). Enriched deuterium (99.96 % fractional ²H abundance; Eurisotop, Cambridge Isotope Laboratories, Saint-Aubin, France) was mixed with the water to obtain a deuterium isotopic composition of ∼500 ‰. To calculate the amount of enriched deuterium to be applied, we used the “D-enriched water” Excel spreadsheet by Alex Sessions of the Division of Geological and Planetary Sciences, Caltech (http://web.gps.caltech.edu/~als/resources, last access: 26 July 2023). To validate the amount of enriched deuterium added, preliminary mixing tests were performed and a sample of the final mixture was sent for isotopic analysis. brilliant blue FCF (Proquimac, Barcelona, Spain) was selected as the dye tracer because of its favourable physicochemical and toxicological characteristics (Flury and Flühler, 1994, 1995). The dye shows a Langmuir-type of sorption behaviour in soil (German-Heins and Flury, 2000). Thus, the higher the concentration, the less relative sorption is expected, as sorption sites become saturated with the dye tracer. In our experiment, a brilliant blue FCF concentration of 5 g L⁻¹ was used (Flury and Flühler, 1994).

2.3 Stemflow infiltration area

Just after the experiment was finished, litter was carefully removed from the surrounding area of the tree base. As per the methodology outlined in Llorens et al. (2022), a plastic mesh was positioned over the soil surface and the infiltration area (m²) was estimated by counting the total number of grid areas covering dye-stained soil and multiplying by the grid-square area (1.27 cm² per grid square).

2.4 Plot excavation and soil sampling

A total of 1 d after the dual-labelled water application, we opened four vertical trenches (30 cm deep and 100 cm wide) in the four cardinal directions at 40 cm distance from the tree bole. Later, two more cross-sections were opened at a distance of 25 and 10 cm in each direction (Fig. 2). Over the course of the day (27 May 2021), a total of 12 soil profiles were excavated and photographed. The face of each profile was levelled by a steel scraper knife to minimize shadows caused by its unevenness. For each of the 12 soil profiles, dye-stained (blue) and non-stained areas were sampled with steel cylinders (100 cm³) to measure the dye concentration and isotopic composition. In total, 63 samples were collected and their location georeferenced. Immediately after the samples were taken, they were covered with plastic film, placed into coolers to maintain soil moisture and prevent evaporation, and transported to the laboratory.

Figure 2 (a) Top view of the experimental plot with squares indicating the location of the vertical soil profiles excavated around the tree. (b) Schematic of vertical soil profiles with dye patterns.

2.5 Image analysis and preferential flow indices

Photographs of the soil profiles were taken using a cell phone equipped with a dual-lens camera system of 12 MP (iPhone 11, Apple Inc., CA, USA). The photographs were edited for geometric correction, contrast and brightness, using Adobe Photoshop 21.2.12 (Adobe Inc., San José, CA, USA). As the photographs were taken at different times of the day with varying illumination (e.g. shadows vs. sunlight), a primary colour scale was utilized to recalibrate the colours and illumination in each photograph (Fig. 3a). The saturation of the blue stains was then maximized to make it easier to distinguish dye intensities (Fig. 3b).

Figure 3 Image analysis procedure of soil profile photographs: (a) geometric correction, recalibration of colours and illumination; (b) colour saturation adjustment; and (c) colour segmentation using a supervised classification.

These corrected images were then used to classify dye-stained and non-stained areas with the maximum-likelihood supervised classification tool in ArcGIS 10.4.1 (Environmental Systems Research Institute Inc., Redlands, CA, USA) (Fig. 3c). The maximum-likelihood classification assigns each cell in the input raster to the class that it has the highest probability of belonging to. Identifiable dye-stained and non-stained areas of interest were manually digitized for each image and used as a training set for the classification tool. Unwanted objects were manually masked out to prevent potential inaccuracies in the classification process.

Digitized images were then further analysed to quantify several preferential flow indices. These indices were calculated as follows:

• Dye coverage (Dc, %) (Flury et al., 1994) is the proportion of the dye-stained area in relation to the total area of the soil profile under consideration. It was calculated in 10 mm depth increments (dye coverage profile) and for the entire profile (total dye coverage). Total dye coverage is lower when there is a higher degree of preferential flow.

• Uniform infiltration depth (UniFr, cm) (Van Schaik, 2009) is the depth at which the dye coverage decreases below 80 %. High values indicate a uniform infiltration process, whereas low values indicate preferential flow.

• Preferential flow fraction (PF_fr, %) (Van Schaik, 2009), defined as the fraction of the total infiltration that flows through preferential flow paths, is calculated as

  $$\begin{array}{}\text{(1)}& {\mathrm{PF}}_{\mathrm{fr}}=\mathrm{100}\cdot \left(\mathrm{1}-{\displaystyle \frac{\mathrm{UniFr}\cdot \mathrm{Wp}}{\mathrm{Tsa}}}\right),\end{array}$$

  where Wp is the width of the profile (cm) and Tsa is the total stained area (cm²).

• Length index (Li, unitless) (Bargués Tobella et al., 2014), defined as the sum of the absolute differences between the dye coverage (Dc, %) of two consecutive depth intervals i and i+1 in a vertical profile of n intervals, is calculated as

  $$\begin{array}{}\text{(2)}& \mathrm{Li}=\sum _{i=\mathrm{1}}^{n-\mathrm{1}}\left|{\mathrm{Dc}}_{i+\mathrm{1}}-{\mathrm{Dc}}_i\right|.\end{array}$$

  This index is related to the degree of heterogeneity of the dye infiltration pattern, with higher Li values indicating greater preferential flow.

• Peak index (Pi, unitless) (Bargués Tobella et al., 2014) is the number of times that the vertical line defined by the total dye coverage intersects the dye coverage profile. This parameter is also related to the heterogeneity of the stained patterns, with high values of the parameter indicating a high degree of preferential flow.

• Mass fractal dimension (D_s, unitless) (Hatano and Booltink, 1992; Hatano et al., 1992) is related to the basic shape of an object (staining patterns), for example a dot, line or plane. D_s can be derived by fractioning a stained profile into squares (pixels) of side r (the smallest value of r corresponding to the width of an original pixel) and then counting the number of partially stained squares, N(r). Repetition of this process with increasing r values gives a series of N(r) values. Double-logarithmic plots of N(r) against r show a linear correlation from which D_s is calculated from the slope as

  $$\begin{array}{}\text{(3)}& \log N\left(r\right)=-D_{\mathrm{s}}\cdot \log r+c,\end{array}$$

  where c is a constant. D_s is calculated in the range of 0 to 2. D_s=0 when the stained pattern is a dot pattern, D_s=1 when it is a line pattern (long and slender) and D_s=2 when it is the dimension of a completely filled area.

2.6 Dye concentration analysis

A similar dye extraction procedure to that of Forrer et al. (2000) was followed. An aliquot of 1 g was taken from the centre of each soil core on the external side. The external side of the soil sampled using the steel ring corresponds to the one captured in the photographs. This aliquot was then mixed with 5 mL of deionized water, shaken overhead for 3 h and centrifuged for 30 min at 3000 rpm (ROTINA 420, Hettich, Kirchlengern, North Rhine-Westphalia, Germany). The concentration of brilliant blue FCF of the supernatant was measured by UV–Vis spectrophotometry at a wavelength of 630 nm, based on a previously established calibration curve relating absorbance (A) to dye concentration (C_d, mg L⁻¹) (C_d=6.83 A; R²=0.99). Because we did not determine the dry mass of the soil cores, the dye concentration in the soil was then expressed as the mass of dye per mass of wet soil (C, mg kg⁻¹ of wet soil).

2.7 Isotopic analysis

A sufficient aliquot to fill a test tube (20 mL) was taken from the internal side of each of the 63 core samples. These samples were sent to the Scientific and Technical Services of the University of Lleida for isotopic analysis. Soil samples were subjected to cryogenic extraction, and the stable water isotopes (δ¹⁸O and δ²H, ‰) were determined by cavity ring-down spectroscopy (Picarro L2120-i analyser, Picarro Inc., Santa Clara, CA, USA). The equipment had an accuracy of <0.1 ‰ for δ¹⁸O and <0.4 ‰ for δ²H, based on the repetition of four reference samples provided by the International Atomic Energy Agency (IAEA). All isotopic data were expressed in terms of δ values and calculated as follows:

$$\begin{array}{}\text{(4)}& \mathit{\delta }=\left({\displaystyle \frac{R_{\mathrm{sample}}}{R_{\mathrm{VSMOW}}}}-\mathrm{1}\right)\cdot \mathrm{1000}\mathrm{‰},\end{array}$$

where VSMOW is Vienna Standard Mean Ocean Water and R is the isotope ratio (¹⁸O${/}^{\mathrm{16}}$O or ²H${/}^{\mathrm{1}}$H).

2.8 Soil characteristics

To characterize the soil at the experimental site, six vertical 30 cm depth profiles were dug in the vicinity (<1 m) around the monitored tree. A depth of 30 cm was the maximum because no significant changes in soil physical properties had been observed in similar forest stands in the same study area beyond that depth (Rubio et al., 2008). Tree roots encountered during excavation were cut and removed. Six soil cores were taken with steel cylinders (100 cm³) in 5 cm depth intervals from the soil surface down to 30 cm in order to determine volumetric soil water content (VWC, cm³ cm⁻³), bulk density (BD, g cm⁻³), soil texture, organic matter content (OM, %) and porosity (ε, %). Once taken, samples were immediately covered with plastic film and placed inside coolers to maintain sample moisture. The volumetric soil water content and bulk density were determined for each core, and the latter was used to calculate porosity, assuming a particle density of 2.65 g cm⁻³ (Hao et al., 2007). The soil organic matter content was determined by loss on ignition, following Wang et al. (2012).

For soil texture analysis, the protocol described in Faé et al. (2019) was used. The dry soil cores were disaggregated and then a representative, well-mixed 50 g aliquot was collected from each sample. Gravel (diameter > 2 mm) was removed by sieving. Organic matter was then removed by use of 15 % hydrogen peroxide, and clay aggregates were dispersed with sodium hexametaphosphate. The soil was next sieved through a 50 µm mesh to separate sand-sized particles (0.05 to 2 mm). To determine the silt and clay fractions, four subsamples were taken from the <50 µm suspension and analysed by laser diffraction (Mastersizer 2000, Malvern Panalytical Ltd., Malvern, Worcestershire, UK). Prior to the analysis, the samples were sonicated for 30 s. Instrument parameters were set to refracting index 1.52, absorption 0.1, stirrer speed 750 rpm, pump speed 1750 rpm and analysis time 30 s. The soil texture (sand, silt and clay fractions) and bulk density were then used to calculate saturated hydraulic conductivity by the Rosetta v.3 pedotransfer functions model (Zhang and Schaap, 2017).

In addition, to determine the isotopic composition of the soil profile before the experiment (see Sect. 2.7), disturbed soil samples were taken in each intercardinal direction from the experimental tree, at a distance of 3 m from the stem. Samples were taken at the following depth increments: 0–5, 5–10, 10–20, 20–30, 40–50 and 90–100 cm.

Figure 4 Soil physical properties at the pine forest stand. The abbreviations and variables used in the figure are as follows: BD – bulk density, K_s – saturated hydraulic conductivity, ε – porosity, OM – organic matter content and δ²H – deuterium isotopic composition of soil water. Note the different scales of the depth axes for the isotopic composition (the green dotted line is the 30 cm depth). Data represent means and standard deviations (n=6 for physical properties and n=4 for isotopic composition).

3 Results and discussion

3.1 Soil characteristics

Data obtained from the 30 cm depth profiles were pooled (n=6) to characterize the physical properties of the soil at the forest stand. Similarly, the isotopic composition values of the four soil cores were also pooled. Figure 4 shows the average and standard deviation of the soil properties. Soil texture was classified as silt loam according to the United States Department of Agriculture (USDA) classification (USDA, 2017). Bulk density increased from about 1 g cm⁻³ at the topsoil to 1.6 g cm⁻³ at 30 cm depth. Similar bulk density values were observed by Rubio et al. (2008) in the same study area at 50 cm depth. We did not expect root growth or development to be restricted at our site, as root growth in silt loam soils will only be affected at bulk density values > 1.6 g cm⁻³ (Zisa et al., 1980). The increase in bulk density corresponds to a decrease in the saturated hydraulic conductivity, porosity and organic matter content with depth (Rubio et al., 2008). Soil porosity decreased by about 24 % from the topsoil to the subsoil, whereas the water content at field capacity declined by only about 2 %. The high organic matter content in the topsoil (16 %) is due to high amounts of litter on the forest floor, as found by Molina et al. (2019). The δ²H decrease (from −48 ‰ in the topsoil to −55 ‰ at 100 cm depth) is due to evaporation, which was also observed by Sprenger et al. (2019) in the same forest stand.

3.2 Stemflow infiltration area

There are few studies that have used dye tracers to calculate the stemflow infiltration area (see the review by Carlyle-Moses et al., 2020). In our experiment, stemflow infiltrated into the soil in an annular shape around the base of the tree trunk and no surface runoff was observed. The area of the annular shape was 0.14 m², which is smaller than the infiltration area of 0.245 m² reported by Schwärzel et al. (2012). However, the other stemflow infiltration areas reported are usually much smaller (e.g. Carlyle-Moses et al., 2018; Metzger et al., 2021; Llorens et al., 2022). In pine forests, there is usually a thick litter layer that repels water (Iovino et al., 2018) and causes the stemflow infiltration area to expand beyond the immediate vicinity of the tree trunk. The stemflow infiltration area found in this study should be interpreted with caution, as it corresponds to a specific amount (50 mm) and discharge (7 L h⁻¹) and to specific antecedent moisture conditions and soil physical properties, and we suspect that it may vary as any of these conditions change.

Figure 5 Time series of the spatial patterns of the volumetric soil water content (SWC) during the experiment. Artificial stemflow was initiated at t=0 min and stopped after around t=390 min (highlighted in blue). Cross marks denote TDR sensors (approximate locations) and the white circle represents the tree trunk.

3.3 Soil water content dynamics

The spatial distribution of the soil water content around the tree before stemflow initiation was highly heterogeneous, and it remained so during the stemflow experiment (Fig. 5). The spatial variability in the SWC showed a decrease in the coefficient of variation throughout the experiment from 17 % to 11 % towards the end. The detailed temporal dynamics of the SWC measured with the different TDR probes is shown in Fig. S1 in the Supplement. A heterogeneous volumetric SWC around a tree, induced by stemflow, was also observed by others (Liang et al., 2011; Tischer et al., 2020). It should be noted that the SWC response caused by stemflow infiltration can change depending on antecedent soil moisture conditions which could lead to different water flow dynamics (Kobayashi and Shimizu, 2007). Table 1 shows the pre- and post-stemflow SWC (0 and 510 min), the maximum SWC peak and the time elapsed until its occurrence for each of the TDR probes.

Table 1 Soil water content measured with TDR probes in different directions around the tree trunk. Pre-stemflow corresponds to t=0 min and post-stemflow to t=510 min. Peak SWC is the maximum during the experiment.

^* No real peak was observed.

The pre-stemflow SWC showed near-saturated conditions for the southeast and west directions at 30 cm distance from the trunk (Figs. 5 and S1, Table 1). These near-saturated conditions explain the lack of response of these TDR probes during the stemflow experiment. The other TDR sensors, however, responded to the stemflow, as shown by the increase in the SWC over time (Figs. 5, S1). As expected, a faster response was observed in the TDR probes closer to the trunk: at 10 cm earlier than at 30 cm (Fig. S1). Nevertheless, the maximum difference in the SWC (i.e. difference between the max SWC and the pre-stemflow SWC) was similar at both distances (Table 1). The post-stemflow SWC (∼2 h after the end of the experiment) showed a small decrease in soil moisture from peak values, indicating that soil water was gradually percolating to deeper soil horizons while the upper horizon remained at field capacity (Fig. S1).

Figure 5 shows that stemflow water was mostly directed to the west, north and east, whereas the southwest and south received less. The northeast and northwest locations at 10 cm distance showed the greatest increases in the SWC (15 % and 16 %), while the SWC in the southeast and south only increased by 5 %. At 30 cm distance, the SWC to the north increased by 19 %, whereas TDR probes to the southwest and west showed no change in the SWC (Figs. 5 and S1, Table 1). In all of the TDR probes that responded to stemflow, the maximum SWC was reached close to the end of the experiment (i.e. around 6 h after initiation of the stemflow).

Our results show that stemflow affects the SWC in the topsoil, even at a horizontal distance of 30 cm from the trunk. Similarly, Carlyle-Moses et al. (2018) used a dye tracer experiment to demonstrate that stemflow infiltrated directly into the soil adjacent to stems and extended belowground to horizontal distances of up to 35 cm from the stems of juvenile lodgepole pine. Contrarily, our finding differs from that of Tischer et al. (2020), who suggested that the topsoil water content at a distance of 20 cm from the trunk of a European beech and sycamore was mainly driven by throughfall rather than stemflow during a stemflow dye tracer experiment conducted on natural rainfall events. On the other hand, Metzger et al. (2017, 2021) found that the SWC was lower in the immediate vicinity of tree stems than in distal areas, due to the differences in soil physical properties (such as lower field capacity and greater saturated hydraulic conductivity near the stem), and suggested that macroporosity enhances drainage near the tree trunk. These dissimilarities may be because our study excluded the throughfall contribution and examined only and exclusively stemflow infiltration dynamics. Our study found that, in the absence of throughfall, the SWC responded similarly to stemflow at 10 and 30 cm from the trunk, indicating that lateral flow was occurring even in the shallow soil layers. We found that stemflow was preferentially funnelled belowground through coarse roots (see Sect. 3.5). Thus, lateral stemflow redistribution within the soil was driven by the root architecture of Scots pines, which develops both horizontal (growing in shallow soil layers) and vertical (growing down to deep soil layers) coarse roots, along with their corresponding fine roots (Fig. S2).

The dynamics of the SWC at different depths also revealed marked differences between the cardinal directions (Fig. 6). In general, the pre-stemflow SWC is greater at depth, and surface probes show a more pronounced response to stemflow towards the end of the experiment (post-stemflow SWC). As expected, the shallower probes exhibited the most pronounced SWC response with a maximum increase ranging from 14 % to 20 %. This is because the surface soil layer (0–10 cm) is the first to become wetter, and the continuous reception of stemflow water contributes to its quicker saturation. However, there is indication of bypass flow, as shown by the faster response of deeper TDR probes, i.e. the TDR probe at 30–60 cm in the east profile and the TDR probe at 30–40 cm in the south profile reacted to the stemflow more quickly than the other probes did (Fig. 6b, c, d). Bypass flow, either by a macropore channel or root presence, leads to spatially irregular wetting of the soil profile with differing flow velocities that allow infiltrating water to move farther and reach deeper soil layers with minimal interaction with the bulk soil (Gerke, 2011).

Figure 6 Temporal evolution of the volumetric soil water content (SWC) in the depth profiles of each cardinal direction during the stemflow experiment. (All times referenced in the figure are local times.)

3.4 Water table dynamics

Water levels monitored by the mini-piezometers indicated the presence of temporary perched water tables around the tree trunk during the stemflow experiment (Fig. 7). Even near-saturated locations (TDR probes swc-SE30 and swc-W30) received stemflow water, as indicated by the occurrence of blue water in nearby mini-piezometers wt-SE20 and wt-W20, respectively (Fig. S3). However, no SWC response to stemflow was observed at these locations because the soil was already saturated. Although some of the mini-piezometers were close to TDR probes, the relationship between them showed some inconsistencies. For example, no water was observed in mini-piezometer wt-NW20, located between TDR probes swc-NW10 and swc-NW30, which showed maximum SWC changes of 15 % and 7 %, respectively, nor was any water found in mini-piezometer wt-N20 situated close to TDR probes swc-N10 and swc-N30, which showed maximum SWC changes of 11 % and 19 %, respectively (Table 1). This result is likely due to water flowing rapidly through the preferential flow channels and not reaching a positive pore water pressure. On the contrary, increasing water levels were observed at mini-piezometer wt-S20, which was located near TDR probes swc-S10 and swc-S30, which both had a maximum SWC change of only 5 %. However, simultaneous increases in water levels and SWC were observed for mini-piezometers wt-NE20 and wt-E20 (to a lesser extent) and TDR probes swc-NE10 and swc-E10, which were close to each other.

Figure 7 Time series of spatial patterns of water level in the 12 mini-piezometers during the experiment. Light pink areas indicate that no water level was detected in the 30 cm deep mini-piezometer, whereas dotted areas indicate zones that were not monitored. Artificial stemflow was initiated at t=0 min and stopped after around t=385 min (highlighted in blue). Cross marks denote mini-piezometers (approximate locations) and the white circle represents the tree trunk.

Mini-piezometers wt-N50, wt-S50 and wt-W50, which were located more than 30 cm from the trunk (Fig. 1), did not register any increase in water level. However, mini-piezometer wt-E50 had a minor increase in water level. Moreover, the groundwater table itself, monitored with northwest and southeast piezometers (1.7 and 2.0 m deep, respectively) at 1.6 m distance from the tree trunk, did not provide any evidence of a change in groundwater levels (Fig. S4). Taken together, the piezometer and SWC measurements down to 30 cm belowground level demonstrate the great temporal and spatial heterogeneity of the soil water dynamics of stemflow infiltration.

3.5 Image analysis of stemflow stain patterns and preferential flow paths

The dye coverage observed in the different vertical soil profiles for the four cardinal directions shows clear indications of preferential flow (Fig. 8). Every cardinal direction shows the same overall pattern of greater dye coverage with increasing proximity to the tree trunk. Similarly, but in horizontal soil profiles, Tischer et al. (2020) reported stemflow stain patterns around the tree base, indicating a decrease in dye coverage with increasing depth from the surface of the organic layer. The distribution of our dye patterns differs from those observed in previous studies of trees conducted on hillslopes (Schwärzel et al., 2012; Spencer and van Meerveld, 2016; Gonzalez-Ollauri et al., 2020), in which dye coverage was found to be lower in the upslope than in the downslope profiles and a clear downward movement of the dye was seen at increasing distance from the tree. We observed the greatest dye coverage in the immediate vicinity of the tree (at 10 cm distance), which suggests that the majority of the stemflow infiltration water was directed vertically. This is likely a consequence of the flat area at our experimental site.

Figure 8 Dye coverage (black denotes soil-stained areas and grey denotes root-stained areas) of the different soil profiles at a distance of (a) 40 cm, (b) 25 cm and (c) 10 cm from the tree trunk. Root-stained areas refer to the area occupied by stained roots in line with the soil profile. The following parameters to quantify the degree of preferential flow are shown within each graph: Dc – dye coverage (%), also shown as a red line; UniFr – uniform infiltration depth (cm); PF_fr – preferential flow fraction (%); Li – length index; Pi – peak index; and D_s – fractal dimension.

The metrics to characterize the different dye patterns (see Sect. 2.5) corroborate the prevalence of preferential flow (Fig. 8, inserts). The maximum dye coverage observed was less than 30 %, highlighting that water flowed through only a small portion of the soil matrix. The metrics of uniform infiltration depth and the preferential flow fraction indicate the absence of matrix flow (UniFr = 0 cm) and the dominance of preferential flow (PF_fr=100 %). The length index (Li) increased with proximity to the trunk and indicated greater irregularity in the dye patterns and, hence, increased preferential flow. Similarly, the peak flow metric (Pi) showed irregular stained profiles. However, lower values were observed at 10 cm distance where direct infiltration is occurring, indicating a less irregular stained pattern. The fractal dimension (D_s) of the stained patterns ranged from 1.79 to 1.96, indicating that the preferential flow pattern occurred in extensive areas (Hatano and Booltink, 1992).

A substantial fraction of bypass and lateral flow occurred along the surface of coarse roots. Living roots can impede flow locally and change the flow towards soil around the root that is non-compacted and highly conductive (Noguchi et al., 1997). In some of the soil profiles, stained roots contributed larger areas to the dye patterns than stained soil (Fig. 8a north, 8b south, and 8c east and west). This pronounced stemflow belowground funnelling along coarse roots is shown in Fig. 9. The same phenomenon was observed in several field studies using dye tracers (e.g. Liang et al., 2011; Schwärzel et al., 2012; Spencer and van Meerveld, 2016; Gonzalez-Ollauri et al., 2020; Tischer et al., 2020) and geophysical techniques (Guo et al., 2020). The dye coverage in our profiles often showed less coverage in the shallow soil layers (0–15 cm) than at greater depth (Fig. 8). This indicates enhanced bypass flow in the shallow soil layers and increasing dye–matrix interaction in deeper layers. This is of great ecohydrological importance, as it may affect processes such as increased moisture, element deposition and contamination in deeper soil layers, where increased dye–matrix interaction has been evidenced. Increased dye–matrix interaction occurs when macropore channels become discontinuous and roots branch out more. This aspect was described by Van Stiphout et al. (1987) as “Internal Catchment”. Moreover, the soil matrix between macropores (i.e. the next smaller pore fraction) must be saturated before water can flow to the next macropore (Beven and Germann, 2013). Although few small stones were found during the excavation, these may influence preferential flow. We observed that some macropores were intersected by stones and extended dye coverage was found around these stones.

Figure 9 Examples of stained roots observed during the excavation of the soil profiles around the tree trunk, illustrating the preferential belowground funnelling of stemflow along coarse roots.

It should be noted that stemflow stain infiltration patterns could behave differently with other tree species and their associated root–soil systems, soil physical and chemical properties, and antecedent soil moisture conditions (Phillips et al., 2019). For example, using dual-tracer infiltration experiments, Luo et al. (2019, 2023) found that the degree of preferential flow was greater in coniferous than deciduous forest soil and that roots enhanced preferential flow, strongly influencing the distribution of the soil water content. This highlights the important role of tree species and forest composition on preferential flows. In addition, slope topography plays a key role in determining stemflow infiltration by diverting most of the stemflow water downslope through gravity, which can lead to overland-flow generation (Gonzalez-Ollauri et al., 2020).

The non-uniform dye patterns and SWC heterogeneity around the tree during the stemflow experiment can be attributed to four factors. First, the preferential flow of stemflow on the trunk itself can affect the aforementioned factors (as shown in Pinos et al., 2021). Although stemflow was evenly discharged around the trunk at a height of 1.3 m, preferential flow on the trunk was observed below this height. Second, the non-uniform infiltration of stemflow into the soil and preferential flow in the soil due to the presence of coarse roots, which redistribute water flow both vertically and horizontally, can cause heterogeneity (Schwärzel et al., 2012). Third, the temporary presence of local perched water tables can result in non-uniform dye patterns and SWC heterogeneity (Liang et al., 2011). Finally, non-uniformity can result from a heterogeneous soil structure (Metzger et al., 2017, 2021).

3.6 Dye and deuterium concentrations in soil samples

As expected, the concentrations of both tracers in the preferential flow paths (i.e. the stained areas) were significantly higher than in the surrounding soil matrix (i.e. the non-stained areas). In general, the concentrations of brilliant blue FCF and of the deuterium isotopic composition in the soil samples decreased with distance from the tree (R²=0.49 for both tracers; Fig. 10). Moreover, the maximum concentrations of both tracers in the different consecutive profiles also decreased with increasing distance from the trunk (Table S1 in the Supplement). The brilliant blue FCF concentrations and deuterium isotopic compositions ranged from 0 to 473 mg kg⁻¹ of soil and from −56.5 ‰ to 277.1 ‰, respectively. For both tracers, the concentrations found in soil were considerably smaller than the original concentration of the labelled stemflow water, with the maximum measured in soil samples representing 9.5 % and 55.4 % of their original concentrations for brilliant blue FCF and deuterium, respectively. This can be explained by (1) sorption of brilliant blue FCF on the trunk bark and soil surfaces and (2) dilution of the stemflow water by residual soil water.

Figure 10 Relationships of the travel distance (three-dimensional coordinate system) with the deuterium isotopic composition and brilliant blue FCF dye concentration for the 63 soil samples. Solid lines indicate logarithmic fits.

Both tracers were significantly correlated (p<0.01, R²=0.78; Fig. S5): there was no evidence that deuterium moved deeper or farther through the soil than the dye tracer. This suggests that the two tracers can be used interchangeably to delineate preferential flow pathways. It is known that brilliant blue FCF sorbs to soil surfaces and is, therefore, retarded to some extent (Flury and Flühler, 1995), but the close correlation of brilliant blue FCF with deuterium in our study shows that the stemflow water was mainly moving along preferential flow pathways with minimal interactions with adsorbing surfaces.

The water samples taken from the mini-piezometer revealed that the temporary perched water tables have a similar brilliant blue FCF concentration to the original stemflow water of ∼5 g L⁻¹ (Table S2). This corroborates the prevalence of preferential flow of the stemflow water once it infiltrates into the soil. However, deeper groundwater (>1 m depth) did not show the presence of dye or deuterium enrichment, indicating that stemflow water did not penetrate deep into the subsoil.

4 Conclusions

Stemflow infiltration occurred in an annulus-shaped area in the vicinity of the trunk and was characterized by high spatiotemporal variability. In some near-stem locations, perched water tables developed at shallow depths. The dye patterns and isotopic composition demonstrated (qualitatively and quantitatively) the dominance of preferential flow in the soil. The belowground funnelling phenomenon was clearly shown by the infiltration of stemflow into the soil along the surface of coarse roots and macropore channels. The root architecture was the main driver of stemflow water redistribution both horizontally and vertically in the soil matrix, followed by bypass flow along the macropore network.

The dual-tracer approach (brilliant blue FCF dye and deuterium) and in situ sensing of soil moisture delineate and quantify stemflow belowground funnelling well. brilliant blue FCF has proven to be an excellent dye tracer to stain and measure preferential flow pathways. The concentration (5 g L⁻¹) used correctly ensured that mobility would not be limited (despite its possible sorption by trunk bark and soil surfaces) by the local conditions of the present study (i.e. rough bark and silt loam soils). On the other hand, deuterium is a conservative tracer, but preferential flow paths can only be identified with extensive soil sampling and analyses.

We conclude that the root architecture of trees leads to preferential flow and, thus, the non-uniform distribution of water in the soil profile. The funnelling effect causes water to infiltrate along roots and macropores near the soil surface and deeper in the soil. Our results are in line with previous research and contribute a more general mechanistic understanding of belowground stemflow infiltration in forests. Most importantly, our results suggests that stemflow belowground funnelling is an important hydrological and ecological feature in humid Mediterranean forest catchments.