Field observations of soil hydrological ﬂow path evolution over 10 Millennia

. The presence or absence of preferential ﬂow strongly controls water ﬂow and transport in soils. It is ubiquitous, but difﬁcult to characterize and predict. This study addresses the occurrence and the evolution of preferential ﬂow during the evolution of landscapes, and here speciﬁcally during the evolution of hillslopes. We targeted a chronosequence of glacial moraines in the Swiss Alps to investigate how water ﬂow paths evolve along with the soil forming processes. Dye tracer irrigation experiments with Brilliant Blue solution (4 g/l) were conducted on four moraines of different ages (30, 160, 3000, 5 and 10000 yrs). At each moraine, three dye tracer experiments were conducted on plots of 1.5 x 1.0 m. The three plots at each moraine were characterized by different vegetation complexities (low, medium, high). Each plot was further divided into three equal subplots for the application of three different irrigation amounts (20, 40, 60 mm) with an average irrigation intensity of 20 mm/h. The day after the experiment ﬁve vertical soil sections were excavated and the stained ﬂow paths were photographed. Digital image analysis was used to derive average inﬁltration depths and ﬂow path characteristics such as the volume and 10 surface density of the dye patterns. Based on the volume density, the observed dye patterns were assigned to speciﬁc ﬂow type categories. The results show a signiﬁcant change in type of preferential ﬂow paths along the chronosequence. The ﬂow types change from a rather homogeneous gravity driven matrix ﬂow in coarse material with high conductivities and a sparse vegetation cover at the youngest moraine to a heterogeneous inﬁltration pattern at the medium-age moraines. Heterogeneous matrix and ﬁnger ﬂow are dominant at these intermediate age classes. At the oldest moraine only macro pore ﬂow It was shown that the complex interaction of vegetation and soil development and its proven effect on ﬂow path development also impacts the water balance, as the storage and conductivity properties of the soil change. However, the interplay between preferential ﬂow paths and the soil hydraulic behavior not only inﬂuences the soil water budget, but also runoff formation. These ﬁndings provide important insights on hydrological ﬂow path evolution in transient systems.

with a reduction in particle size and an increase in porosity leading to an increase in subsurface water storage. To test these hypotheses we targeted a chronosequence of glacial recession in the Swiss Alps.
Areas with receding glaciers have been shown to be suitable for soil development studies [ Crocker and Dickson (1957), Douglass and Bockheim (2006), He and Tang (2008), Egli et al. (2010), Dümig et al. (2011), Vilmundardóttir et al. (2014), D'Amico et al. (2014)]. In the cool and humid climate regions of former glacial areas the soils develop from mineral soils 15 to soils with a highly organic topsoil. These organic soil types are less intensively studied with regard to their soil hydraulic behavior compared to mineral soils (Carey et al., 2007). It is known that these soils differ in their soil hydraulic properties from mineral soils (high total porosity (up to 90%) and a low bulk density (Carey et al., 2007)) but little is known about how this development impacts water flow paths. Dye tracer experiments, and an analysis of soil texture and soil physical properties were used to investigate how water flow paths evolve with hillslope age. Dye tracer experiments combined with digital image 20 processing have been applied successfully to study preferential infiltration in soils [Weiler (2001), Bogner et al. (2008), Blume et al. (2008, Laine-Kaulio et al. (2015), Hardie et al. (2011), Cheng et al. (2014)]. In our study we use this method to identify how flow paths change during the co-evolution of soil, vegetation and topography. Understanding the changes in preferential flow paths as a result of the natural co-evolution of landscape forming factors can provide valuable knowledge on how these systems can also change as a result of human intervention (Richter and Mobley (2009)). 25 2 Material and methods

Study site
The study area is located in the foreland of the Stein glacier above the tree line in the Central Swiss Alps, south of the Sustenpass in the Urner Alps (appr. 47 • 43'N, 8 • 25'E). Elevations range from 1900-2100 m a.s.l. The area lies in the polymetamorphic "Erstfelder" gneiss-zone, which is part of the Aar-massif (Blass et al., 2003). The geology is defined by metamorphosed pre- 30 Mesozoic, metagranitoids, gneisses, and amphibolites [Heikkinen and Fogelberg (1980), Schimmelpfennig et al. (2014)], thus the material is mainly acidic and rich in silicate. The closest official weather station is located 18 km away at Grimsel Hospiz (46 • 34'N, 8 • 19'E) at an elevation of 1980 m a.s.l. For the norm period from 1981 to 2010 the station recorded an annual mean temperature of 1.9 • C and an annual precipitation of 1856 mm. The precipitation distribution throughout the year is fairly uniform with a slight increase in the winter months (Schweizerische Eidgenossenschaft, 2016). The glacier foreland consists of moraines with unconsolidated glacial till. The humid and cool climate together with the nutrient-poor substrate and a relative high water permeability of the glacial till favor the formation of podsolic soils and humus in this area (Heikkinen and Fogelberg, 1980). 5 The moraines of the Stein glacier were exposed due to its retreat to the south. Four moraines were selected for this study (see Fig. 1). Schimmelpfennig et al. (2014) conducted a detailed dating study of the Stein glacier moraines, based on high-sensitivity beryllium-10 moraine dating and found that the ages of three moraines range between 160 to 10,000 years. The age of a fourth moraine was dated to 30 years based on maps and aerial photos.  The oldest moraine with an age of 10,000 years is facing north-east. The second oldest moraine with an age between 2,000 10 and 3,000 years is the only one facing south. The two youngest moraines were exposed in the years 1860 and 1980-1990, and thus have an age of 160 and 30 years, respectively. Both moraines are facing north-east and are located closer to the glacier tongue at a distance of approximately 1 km from the oldest moraines. Both moraines are south of the glacial lake "Steinsee" (1930 m a.s.l.), which is a proglacial lake that was formed by the glacier retreat in 1924 (Blass et al., 2003).
The vegetation of the moraines was mapped in summer 2017 by the project partners of the Geobotany Group of the University 15 of Freiburg, Germany (Maier et al., 2019). Pronounced differences in vegetation coverage and species distributions were found among the four age classes. The vegetation of the oldest moraine was dominated by a variety of prostrate shrubs together with small trees and several grasses. On the 3,000-year-old moraine a grassland cover with fern, mosses, sedges and forbs was found. The two youngest moraines, however showed a lower degree of vegetation complexity. On the 160-year-old moraine a combination of grasses, lichen, forbs, and shrubs was present. The youngest moraine still shows only a sparse vegetation cover 20 4 https://doi.org/10.5194/hess-2020-28 Preprint. Discussion started: 24 January 2020 c Author(s) 2020. CC BY 4.0 License. with mainly grass, moss, forbs, and a few shrubs.
Three study plots were selected at each moraine, based on degree of vegetation complexity (low, medium and high complexity).
Vegetation complexity is characterized by the plant structural complexity, which is determined by plant species diversity. A more detailed description of the application of this method can be found in Musso et al. (2019).

Field experiments 5
The dye tracer experiments were conducted between mid July and mid August 2018. We used Brilliant Blue FCF as dye tracer due to its good visibility, high mobility and non-toxicity. We used a concentration of 4 g/l at which the tracer shows good sorption and visibility (Weiler and Flühler, 2004). The size of each study plot was 1.5 x 1.0 m. Each plot was further divided into three equal subplots of 0.5 x 1.0 m. The subplots were irrigated with three different amounts of dyed water (20 mm, 40 mm, 60 mm) and an irrigation intensity of 20 mm/h (Fig. 2). The irrigation intensity is relatively high with a return period of 2.8 years (Fukutome et al., 2017). In preparation of the tracer application large vegetation in form of shrubs and bushes were cut off to a height of a few centimeters to reduce interception. The tracer was applied with a hand-operated sprayer connected to a battery powered pump which guaranteed a constant pressure for a uniform flow rate of 60 l/h. The tracer was applied in three steps. In each step, the alternating intervals of irrigation and breaks varied. In the first step all three subplots were irrigated for 60 minutes in a sequence of 5 minutes 15 irrigation and 5 minutes break. After finishing the first step the first subplot was covered for protection. In the second step, the other two subplots were irrigated for additional 60 min in the same sequence of 5 min irrigation and a 10 min break. In the last step only the third subplot was irrigated for 60 min in a sequence of 2 min irrigation and 10 min breaks while the other two plots remained covered. After the tracer application the whole plot was covered. The next day the plot was excavated in up to five profiles of 7 to 10 centimeters. After the profile cuts were made with pickaxes, spades and hand shovels, the profile walls were cleaned. Hanging roots were cut off and rocks were not removed but made visible. The photographs of the profile walls were taken with a Panasonic Lumix DMC-FZ18 camera and a resolution of 2248 x 3264 pixels. A big umbrella was used to provide a uniform light distribution in the photographs and to avoid direct sunlight. A wooden frame for a geometric correction 5 and a gray-scale (Kodak) attached to the frame (Fig. 2) for a later color adjustment were included in the photographs.

Image analysis
The image analysis procedure by Weiler (2001) was used to generate binary images of the photographs showing stained and unstained areas (Fig. 2). A detailed description of the method can be found in Weiler and Flühler (2004). Instead of the original IDL software package a similar Python version was used. Basically, a geometric correction, a background subtraction and color 10 adjustment was carried out to correct differences in image illumination and changes in the spectral composition of daylight.
The delineation of rocks and plants was done manually. In the resulting binary image the horizontal and vertical length of a pixel correspond to 1 mm.
For a quantitative comparison of the dye patterns the maximum infiltration depth, the volume density and surface area density, as well as the stained path width were calculated. These parameter are frequently used for an objective comparison and surface area density are originally steorological parameters which are used to relate three-dimensional structures to measured two-dimensional parameters (Weibel, 1979). The volume density corresponds to the dye coverage and can be derived from one-dimensional information by calculating the fraction of stained pixels for each depth. The volume density profile is defined 20 by the fraction of stained pixels per depth and is calculated as the average of all excavated profiles per plot. The surface area density in one-dimension is calculated by using the intercept density, which describes the amount of intercepts between stained and unstained pixels divided by the horizontal width of the soil profile. The profile of the surface area density describes the amount of intercepts per depth and is then also averaged over all photographed profiles per plot. Volume density provides no information whether the stained area is the sum of many small fragments or a few large ones, thus the volume density alone 25 should not be used to characterize flow patterns and the surface area density should be used as a supplementary parameter. A high surface area density indicates a large number of small features.
Following the method described by Weiler (2001) the resulting dye patterns were next classified into flow type categories based on the proportions of three selected stained path width classes (stained path width <20 mm, 20 mm-200 mm, >200 mm) relative to the volume density. The stained path width is equal to the horizontal extent of a stained flow path (Weiler, 2001). 30 This classification method distinguishes between five flow types: (1) macro pore flow with low interaction, (2) mixed macro pore flow (low and high interaction), (3) macro pore flow with high interaction, (4) heterogeneous matrix flow/finger flow, and (5) homogeneous matrix flow. Dye patterns, which cannot be classified as one of these flow types are categorized as undefined.
This method is based on the assumption that the dye patterns are mainly controlled by certain preferential flow processes and each flow process creates a characteristic dye pattern that can be described by the extent and distribution of stained features.
This method was proven to be suitable for the investigation of Weiler (2001) and Weiler and Flühler (2004) and was also used in a variety of additional studies [Bachmair et al. (2009), Gimbel et al. (2016, Mooney and Morris (2008)]. However, an extension of the flow type categorization was needed for our study site with soils of different ages, texture and additionally a high stone content. In the extended classification, we avoid a clear differentiation between macro pore flow and finger flow.

5
The original classification assigns finger flow only when both, the medium sized stained path width class (20-200 mm) and the class with the biggest stained path widths (> 200 mm) account for approximately half of the dye coverage. This implies that finger flow is only prevalent when the dye pattern is characterized by a majority of broader stained path widths, ignoring that the size of the finger-like flow paths can vary over a broad range (Wang et al., 2018a). Thus, we argue that while finger flow and macro pore flow are caused by different properties, they can lead to similar dye patterns and distributions of stained path 10 width classes. This is especially the case for macro pore flow with high interaction, which creates broader stained paths that could also be assigned to finger flow. Therefore, both flow types were considered in the extended classification for this class.
Furthermore, it was observed that the presence of rocks within the image analysis interrupts homogeneous blue stained areas and thus leads to smaller stained path widths. Using the original classification scheme on a soil profile with high stone content suggests a heterogeneous flow pattern, which can be classified as heterogeneous matrix flow, finger flow, or macro pore flow 15 depending on the abundance of rocks. Therefore, an additional class has been introduced, which is used when homogeneous matrix flow between rocks takes place. The classification rule for the additional flow type class is based on the proportion of blue dye coverage of the available permeable matrix space (profile width minus sum of stone widths per row). If at least 95 % of the permeable space is stained by blue dye the flow type is classified as matrix flow between rocks.

Soil sampling and laboratory analysis 20
Soil samples were taken during August and September of 2018 close to each dye tracer plot. For particle size analysis, two disturbed soil samples per depth were taken at 10, 30, and 50 cm depth at each plot. The total of 72 samples was analyzed in the laboratory between November 2018 and January 2019 by using a combination of dry sieving (particles > 0.063 mm) and sedimentation analysis (particles < 0.063 mm) with the hydrometer method (Casagrande, 1934). Organic matter removal was only possible by floating off the lighter fractions prior to particle size analysis. Since three plots were selected per moraine for 25 the brilliant blue experiments, six samples per depth and age class (a total of 18 samples for each moraine) were available for the particle size analysis. Particles between 2 mm and 0.063 mm were classified as sand, between 0.063 mm and 0.002 mm as silt and particles smaller than 0.002 mm as clay.
Furthermore, at each plot two 250 cm 3 and one 100 cm 3 undisturbed soil samples were taken at a depth of 10 and 30 cm.
Three samples of 100 cm 3 were taken at 50 cm depth. Thus, per age class nine samples per depth were available for the 30 determination of the structural parameters porosity and bulk density. The porosity was determined in the lab using the water saturation method.

Statistical analysis
The non-parametric Kruskal-Wallis test was used to test the significance of the differences in the soil texture among the four age classes. It can be applied when the assumption of a normal distribution can not be made and is also valid for small sample sizes. We applied the test to each particle fraction across the four age classes. Average values were based on 18 samples per age class. The grain size distribution at 10 cm depth of the oldest moraine was excluded from consideration, as due to the high 5 organic content not the entire organic matter could be removed and the results may therefore be erroneous.

Soil texture and structural parameters
Comparing the depth averaged soil texture over the millennia we find that while the soil texture at the youngest moraine mainly consists of sand, the particle sizes decrease over the millennia with silt being the largest fraction after 10,000 years ( Figure 3). 10 Clay content increased with age for the three older moraines, with the youngest moraine being the exception (having the second highest clay content). The Kruskal-Wallis test with a 0.05 confidence level showed that differences in particle size fractions among the four age classes were statistically significant (p-values < 0.05; sand: p=0.0013, silt: p=0.0006, clay: p=0.0018).     The porosity observed at the youngest moraine ranges between 0.22 and 0.37, with no pronounced differences among the three soil depths. The 160-year-old moraine has a higher porosity in the upper 10 cm than the 30-year-old moraine. After 3,000 years the increase in porosity continues but is now also visible in 30 and 50 cm, but with much higher values at 10 cm ( Figure   4a). After 10,000 years the porosity at 10 cm ranges from 0.6 to up to more than 0.8. The other two depths also experienced a further increase in porosity. The decrease in bulk density is also most pronounced in the top layer of the soil (Figure 4b). While 5 after 30 years the bulk density in the upper 10 cm ranges around 1.7 g/cm 3 the bulk density after 10,000 years is much smaller and ranges between 0.2 and 0.7 g/cm 3 . After 3.000 years the trend is also visible in 30 and 50 cm. 9 https://doi.org/10.5194/hess-2020-28 Preprint. Discussion started: 24 January 2020 c Author(s) 2020. CC BY 4.0 License.
Even though the differences in soil physical characteristics between 30 and 160 years are comparatively small, the rates of change during this initial phase are highest (Table 1). Between 160 and 3,000 years, the rates of change are significantly reduced and remain in a similar range between 3,000 and 10,000 years.

Vertical dye pattern analysis
Flow patterns traced with brilliant blue dye changed considerably with hillslope age ( Figure 5). The volume density profile is a 5 measure for the amount of blue dye per depth. The profile patterns of volume and surface area density show distinct differences among age groups, while differences between the vegetation complexity levels are not as clear ( Figure 5 and Figure 6).   The volume density of the blue dye was classified in three selected stained path width groups (Weiler, 2001). Additionally, the volume density of rocks was also classified in three groups ( Figure 5).
An analysis of the average or maximum infiltration depth based on the dye profiles was not possible because not all profiles 11 https://doi.org/10.5194/hess-2020-28 Preprint. Discussion started: 24 January 2020 c Author(s) 2020. CC BY 4.0 License.
could be excavated up to the maximum infiltration depth. In most cases large boulders prevented further excavation or the infiltration depth was more than one meter. The latter was mostly the case at the youngest moraines. Only at the oldest moraine a maximum infiltration depth could be determined based on the dye profiles.  Comparing the 160-year-old moraine to the youngest moraine we find that the volume density is lower and the surface density is higher in the upper part of the profile. Also unstained areas (colored white in Figure 5) are visible, which indicates that the higher surface density is not caused by the existence of rocks that split up a homogeneously stained area as it was the case for the youngest moraine. In this case there is no total dye coverage of the permeable soil and the preferential flow paths are initialized already near/at the soil surface. The surface density profiles show a decrease in the lower half of the soil profile, 15 which either goes along with a decrease in blue dye coverage and an increase in stone coverage or an increase in blue dye coverage without a significant change in stone coverage. Both indicates a reduction in the amount of separate flow paths and an increase in flow path widths, even if the permeable space is reduced due to an increase in rock content. Also the fraction of stained path widths (SPW) bigger than 200 mm increases, which indicates that the dye plumes widens in deeper soil depths. 20 Compared to the two youngest moraines the moraine of 3,000 years shows in general a higher surface area density and a lower dye coverage combined with a higher fraction of unstained permeable soil matrix. This indicates that similar to the 160 years old moraine water is transported in individual flow paths, but here there are more flow paths and they have a smaller width. This can also be seen in the less frequent appearance of stained path width higher than 200 mm ( Figure 5). Similarly to the 160-year-old moraine, preferential flow paths are already initialized at the top of the soil during infiltration (apparent from 25 the white colored areas across the profile in Figure 5).
The oldest moraine with an age of 10,000 years shows the highest surface area density and the lowest volume density combined with the lowest infiltration depths. The surface and volume density profiles show the same pattern: After a peak close to the soil surface, both density profiles show a decrease with soil depth. This means that the dyed water is only transported deeper into the soil via a few individual flow paths. Using the information in the volume density profiles and the stained path widths to characterize flow types (Weiler, 2001) we found a trend from a rather homogeneous flow pattern with gravity driven matrix flow, especially in the top soil, at the youngest moraine to a more heterogeneous flow pattern with a mix of heterogeneous matrix flow and finger flow at both medium age moraines (Figure 7). While the flow characteristics of both the 160 and 3,000-year-old moraine are dominated by finger flow with smaller stained path widths, this is much more pronounced in the 3000-year-old moraine. By contrast, the oldest moraine 5 stands out clearly from the other age groups. Here, only macro pore flow is predominant. For the deeper soil layers, this was confirmed during the field experiments, since the few macro pore flow paths were clearly visible. For the topsoil, the result is less certain, as the blue areas were very difficult to identify during the image analysis due to the very dark color of the organic layer.
The relative frequency distribution of the flow types per moraine age class derived from the results in Figure 7 show a clear 10 shift in the flow type distribution along the age classes (Figure 8). At the youngest moraine all types of finger flow and matrix flow are present and the frequency distribution does not show a distinct peak at any flow type. With increasing age macropore flow becomes more and more important and the peaks in the frequency distribution become more and more pronounced ( Figure   8).

Evolution of soil texture and structure
The observation of bulk density, porosity and soil texture show significant differences between the age groups as well as some clear trends with age. The 30-year-old soil is characterized by coarse material, with a soil texture composed of almost three quarters of sand. The soil texture observed at the 160-year-old moraine does not differ strongly from the 30-year old moraine, 5 whereas slight changes are already evident in porosity and bulk density. These findings are similar to findings by Dümig et al. (2011), who found no specific trend in particle size distribution for soils in the range of 15 to 140 years at a soil chronosequence in the foreland of the retreating Damma glacier (Switzerland). They furthermore found a high variability in particle size distribution within the single age classes. However, in the same study a slight decrease in variability with increasing age and a noticeable higher clay content was found at the reference site with an age of more than 700 years. He and Tang (2008) also 10 revealed a non-linear increase in maximum clay content for soils up to 180 years at a glacier foreland in a monsoon-temperate region in southwest China.
At 3,000 years of soil development we observe a distinct increase in silt and a reduction in sand content and this development continues, as observed in the 10,000 year moraine, where the silt content now makes up the largest share. These findings correspond with findings by Douglass and Bockheim (2006) who studied several moraines in Buenos Aires with ages ranging from 15 16,000 year to even 1,000,000 years and found an accumulation of clay-sized particles with increasing age, but with a decrease in accumulation rate over the years.
The soil material at the 30-years-old moraine showed a relatively uniform porosity and bulk density throughout the profile.
After 160 years of soil development the porosity in the top layer increased. As these changes cannot be traced back to changes in the grain size composition, the increase in porosity is probably due to the development of a more dense vegetation. The

25
The bulk density shows a reduction in the top layer of the 160-year-old moraine. The decrease in bulk density is also linked to the change of particle sizes and the development of vegetation with time, which includes accumulation of organic matter and the development of a root system. Similar findings in bulk density evolution were also observed by Crocker and Major (1955) who found a decrease in bulk density over the first 200 years of soil development from more than 1.4 g/cm 3 to less than 0.8 g/cm 3 for glacial till in southeastern Alaska. A less pronounced reduction was also found by Crocker and Dickson (1957). He 30 and Tang (2008) found a reduction for the time span of 180 years from appr. 1.42 to 0.95 g/cm 3 that was also more distinct in the upper horizon. Vilmundardóttir et al. (2014) revealed at a glacier foreland in southeast Iceland under maritime climate conditions a reduction from 1.36 to 1.07 g/cm 3 for a time span of 120 years . All studies mentioned above linked this decrease in bulk density to the vegetation development with time.
After 3,000 years, porosity increases and bulk density decreased even further, and this development is now also visible in deeper soil depths. The continuous increase in porosity and reduction in bulk density can still be attributed to the continuing change in soil texture on the one hand and on the other hand to the pronounced vegetation development and the resulting accumulation of soil organic matter with an even more dense root network that is now over 35 cm deep.

5
The oldest moraine shows a significantly higher porosity compared to the 3,000-year-old moraine and a significantly lower bulk density. The change is visible at all soil depths, with the porosity in the uppermost depth being distinctly higher than the other depths. These differences in soil properties between the soil layers also indicate a progressive formation of distinct horizons in the soil. The significantly higher porosity in the upper layer of the oldest moraine is caused by its thick organic layer (thickness up to 20 cm), which is characterized by porosity of up to 90 percent [this was also found by Nyberg (1995) 10 in sandy-silty till on the west coast of Sweden and Carey et al. (2007) in organic soils in a permafrost region in northwest Canada]. Musso et al. (2019) investigated the evolution of pore sizes in the top 5 cm at the same soil chronosequence and found an increase in number of small soil pores and a decrease in relative proportion of macropores (pore diameter > 0.05 mm) between 160 and 10,000 years. Thus the high porosity in the organic top layer at the oldest moraine is mainly composed of small pores.

15
The top layer therefore has an increased water storage and water holding capacity. Due to the finer soil texture and higher porosity the total storage water capacity of the oldest moraine is larger than that of the the younger moraines. An investigation of the saturated hydraulic conductivity evolution of the near-surface (in 0-5, 5-20, and 20-40 cm) at the same chronosequence by Maier et al. (2019) found a decrease with increasing moraine age and soil depth. The reduction in saturated conductivity was found to be positively correlated with soil texture, indicating that the increasing fraction of fine particles with increasing 20 age have even a bigger effect on the saturated conductivity evolution than the root network development (Maier et al., 2019).

Evolution of flow paths
The flow type classification by Weiler (2001) was used to classify the volume density patterns into flow type categories. A comparison between the observations made during the excavation and the derived flow types showed that in this case an adaptation of the flow type classification was necessary. On the one hand, this adaptation involved the treatment of rocks to 25 prevent mis-classification and, on the other hand, we introduced the possibility of finger-like flow paths with smaller widths.
After these adaptions the derived flow types correspond well with the observations made in the field.
The observed staining patterns and derived flow types show a significant difference between the age groups, whereas no significant difference was observed with respect to vegetation complexity and irrigation amount. At the 30-year-old moraine the water infiltrates homogeneously into the soil, probably due to the very low vegetation coverage. The dye pattern showed a 30 mainly homogeneous staining of the soil material, thus derived flow types are mainly matrix flow in form of homogeneous and heterogeneous matrix flow as well as matrix flow between rocks. Also finger flow occurs at the boundary to the clay layer or is caused by large blocks of rock, which are surrounded by clay. The determined macro pore flow takes place only within the clay layer at a depth below 50 cm. In the clay layer, no significant macropores were identifiable, which is why it is assumed that the water is transported in cracks or along material interfaces. The upper coarse soil material with large pores and a low water holding capacity causes the water to be transported quickly deeper into the soil. The water flow is mainly driven by gravity.
After 160 years the derived predominant flow types shift to heterogeneous matrix flow and finger flow. The observed widening of the dye plumes in deeper soil depths might be caused by a change in material or a reduction of hydrophobicity with soil depth where the influence of plants and organic material decreases (Blume et al., 2009). The dye coverage images show unstained 5 soil areas starting also at the top of the soil profile, which indicates that preferential flow paths are initialized already at the soil surface or in the near-surface layer. This was also observed during the irrigation where we saw that the irrigated water often mainly infiltrated in depressions. Grass patches also tend to inhibit infiltration. It was observed that in the presence of dense grass patches, the water infiltrated only next to the patches, leaving the area below the patches unstained.
Similar observations were made at the 3,000-year-old moraine. The image analysis revealed that preferential flow paths start 10 at the soil surface and thus are initiated by vegetation and micro topography causing a heterogeneous infiltration pattern. Field observations also revealed that heterogeneous infiltration was not only created due to dense grass patches, but also occured under relatively homogeneous grass cover. A laboratory test of the water drop penetration time [DeBano (1981), Doerr et al. (2000)] on a soil sample of the upper soil material showed that the organic layer is highly water repellent in dry conditions (air dried for 2 weeks, water drop penetration time > 10 minutes). An increase in the Hydrophobicity Index (Tillman et al., 1989) 15 with increasing moraine was also found by Maier et al. (2019). Thus, we conclude that hydrophobicity of the organic top layer has a big impact on infiltration and the initiation of unstable flow. Compared to the 160-year-old moraine the 3,000-year-old moraine is characterized by a higher number of narrower preferential flow paths.
The derived dominant flow type class at the 3,000-year-old moraine is macro pore flow with high interaction/finger flow. Of both possible flow types, finger flow is the prevalent flow process causing the dye pattern. Several studies linked the formation 20 of finger-like flow paths to hydrophobic properties of the soil [Wallach and Jortzick (2008), Dekker and Ritsema (2000), Ritsema and Dekker (1994), Blume et al. (2008), Wang et al. (2018a), Hardie et al. (2011)]. It is assumed that hydrophobic compounds that are released during the decay of litter (Reeder and Jurgensen, 1979) or by root activity (Doerr et al., 1998) coat soil particles or are deposited in the pore space and thus create a hydrophobic soil matrix (Doerr et al., 2000). The humid and cool climate of former glacial areas leads to a slow decomposition of vegetation and thus to an accumulation of hydrophobic 25 compounds (Doerr et al., 2000).
At the oldest moraine, we see a significantly lower infiltration depth ( Figure 5). During the experiment no surface runoff was observed. Most of the water was stored in the organic top layer. The soil beneath the top layer was almost completely unstained and water was transported only via a few macropores into deeper layers. A dense network of roots was only observed in the organic top layer, which included the thicker roots of the alpenrose (Rhododendron ferrugineum). The root network in the soil 30 underneath was less dense with roots of smaller diameters but extended to a depth of more than 50 cm. Although the vegetation cover has been reduced to lower interception, the interception storage capacity at the oldest moraine is still comparatively high.
Thus, a reduction in the water available for infiltration cannot be ruled out.
The rock content at the 30, 160, and 3,000-year-old moraines is relatively high, with especially large rocks (widths > 20 mm) in deeper parts of the soil ( Figure 5). The large rocks lead to a reduction of the permeable area and thus can cause funnel 35 flow (Hendrickx and Flury, 2001). This type of preferential flow was especially observed at the youngest moraine, where large boulders (> 25 cm) located at deeper soil depths were surrounded by unstained fine textured material. Smaller sized rocks in the upper part seemed not to have an influence on water transport (apart from reducing the flow-through volume), since these rocks and the surrounding soil were completely stained.
A splitting of flow paths caused by rocks was also observed a few times at the 160-and 3,000-year old moraines. In this case, 5 water flowing past the sides of medium to large sized rocks create a type of finger flow that is not caused by water repellency or air entrapment. Integrating all of our findings on soil structural parameters, texture, vegetation cover and flow path patterns provides an overview over their co-evolution and highlights the derived major flow path controls (Figure 9).  Figure 9. Sequence of observed characteristic dye patterns and derived flow path controls compared to the qualitative evolution of soil texture, structural parameters and vegetation. The shade of the root mass-distribution-triangles is a measure for the vertical root mass distribution with darker color indicating a higher root mass. The width of the triangles is a measure for the root mass comparison between the moraine ages, with broader triangles indicating a higher root mass.
Along the co-evolution of soil and vegetation over 10,000 years the major controls of subsurface flow paths change. At the youngest moraine flow paths are only controlled by soil texture. The coarse material leads to mainly gravitation driven flow.
Preferential flow paths only occur at the interfaces between coarse and finer material.
At the medium aged moraines flow paths are mainly controlled by vegetation shielding, microtopography and hydrophobicity.
The latter is assumed to have an increased impact at the 3,000-year-old moraine. After 10,000 years of hillslope evolution 5 subsurface water transport is highly preferential and controlled by flow paths caused by root channels or boundaries of textural classes. Water storage in the organic layer which is also the main rooting horizon increases strongly.

Conclusions
Using brilliant blue dye experiments and soil sampling we investigated the evolution of water transport paths along soil forming processes. To our knowledge this is the first study examining flow path evolution across the millenia in such detail. The 10 evolution of the grain size distribution shows that particle size decreases with increasing age. The biggest changes are in the sand and silt fraction. Furthermore, also water flow defining structural parameters such as porosity and bulk density change during soil development, resulting in an increasing water storage capacity with age. The depth dependent evolution of these parameters confirms our hypothesis that the soil material develops with increasing age from a homogeneous structure to a layered soil with vertical gradients in flow and storage defining soil properties. Changes in these flow defining parameters are 15 caused by the evolution of grain size distribution and vegetation.
We also confirmed the hypothesis that subsurface flow paths change significantly through the millenia, which is indicated by the derived flow types. Flow types change from homogeneous gravity driven matrix flow to a heterogeneous matrix and finger flow over the first 100-1,000 years. At very young moraines the water is homogeneously distributed within the soil matrix.
However, the water storage capacity is relatively low due to the coarse material and water is transported quickly deeper into 20 the soil due to gravity driven flow. At the medium age moraines, water is transported preferentially via finger-like flow paths deeper into the soil by leaving parts of the soil dry.
With increasing hillslope age, we expected macropores induced by root activities to become more important. After 10,000 years, were the amount of soil matrix macropores decreased significantly, macro pore flow along roots plays an important role, but is not very pronounced. Only a few roots reach beyond the organic top layer. However, this allows a fast transport of water 25 from the upper layer into deeper soil. The organic top layer has a pronounced influence on the soil water budget, by storing a significant amount of water. The increase in water storage with increasing age of the moraines also caused a reduction in infiltration depth.
The influence of preferential flow paths increases with soil age. Preferential flow is, however, not only caused by macropores, but especially for the medium-age moraines seems more controlled by soil surface characteristics such as vegetation patches,  Blume, T., Zehe, E., and Bronstert, A.: Investigation of runoff generation in a pristine, poorly gauged catchment in the Chilean Andes