Characterization of stony soils ’ hydraulic conductivity using 1 laboratory and numerical experiments 2

10 Determining soil hydraulic properties is of major concern in various fields of study. Although 11 stony soils are widespread across the globe, most studies deal with gravel-free soils so that the 12 literature describing the impact of stones on the hydraulic conductivity of a soil is still rather 13 scarce. Most frequently, models characterizing the saturated hydraulic conductivity of stony 14 soils assume that the only effect of rock fragments is to reduce the volume available for water 15 flow and therefore they predict a decrease in hydraulic conductivity with an increasing 16 stoniness. The objective of this study is to assess the effect of rock fragments on the saturated 17 and unsaturated hydraulic conductivity. This was done by means of laboratory experiments and 18 numerical simulations involving different amounts and types of coarse fragments. We 19 compared our results with values predicted by the aforementioned predictive models. Our study 20 suggests that considering that stones only reduce the volume available for water flow might be 21 ill-founded. We pointed out several drivers of the saturated hydraulic conductivity of stony 22 soils, not considered by these models. On the one hand, the shape and the size of inclusions 23 may substantially affect the hydraulic conductivity. On the other hand, laboratory experiments 24 show that an increasing stone content can counteract and even overcome the effect of a reduced 25 volume in some cases: we observed an increase in saturated hydraulic conductivity with 26 volume of inclusions. These differences are mainly important near to saturation. However, 27 comparison of results from predictive models and our experiments in unsaturated conditions 28 shows that models and data agree on a decrease in hydraulic conductivity with stone content, 29

soils, not considered by these models. On the one hand, the shape and the size of inclusions 23 may substantially affect the hydraulic conductivity. On the other hand, laboratory experiments 24 show that an increasing stone content can counteract and even overcome the effect of a reduced 25 volume in some cases: we observed an increase in saturated hydraulic conductivity with 26 volume of inclusions. These differences are mainly important near to saturation. However, 27 comparison of results from predictive models and our experiments in unsaturated conditions 28 shows that models and data agree on a decrease in hydraulic conductivity with stone content, 29

Introduction 1
Determining soil hydraulic properties is of primary importance in various fields of study such 2 as soil physics, hydrology, ecology and agronomy. Information on hydraulic properties is 3 essential to model infiltration and runoff, to quantify groundwater recharge, to simulate the 4 movement of water and pollutants in the vadose zone, etc. (Bouwer and Rice, 1984). Most 5 unsaturated flow studies characterize the hydraulic properties of the fine fraction (particles 6 smaller than 2 mm of diameter) of supposedly uniform soils only (Bouwer and Rice, 1984;7 Buchter et al., 1994;Gusev and Novák, 2007). Nevertheless, in reality, soils are heterogeneous 8 media and may contain coarse inclusions (stones) of various sizes and shapes. 9 Stony soils are widespread across the globe (Ma and Shao, 2008) and represent a significant 10 part of the agricultural land (Miller and Guthrie, 1984). Furthermore, their usage tends to 11 increase because of erosion and cultivation of marginal lands (García-Ruiz, 2010). Yet little 12 attention has been paid to the effects of the coarser fraction on soil hydraulic characteristics, so 13 that the relevant literature is still rather scarce (Ma and Shao, 2008;Novák and Šurda, 2010;14 Poesen and Lavee, 1994). 15 Many authors consider that the reduction of volume available for water flow is the only effect 16 of stones on hydraulic conductivity. This hypothesis has led to models linking the hydraulic 17 conductivity of the fine earth to those of the stony soils. They predict a decrease in saturated 18 hydraulic conductivity of stony soil ( K se ) with an increasing volumetric stoniness ( R v ) 19 (Bouwer and Rice, 1984;Brakensiek et al., 1986;Corring and Churchill, 1961;Hlaváčiková 20 and Novák, 2014;Peck and Watson, 1979; Ravina and Magier, 1984). 21 However, a number of studies do not observe this simple relationship between the hydraulic 22 conductivity and the stoniness (Zhou et al., 2009;Ma et al., 2010;Russo, 1983;Sauer and 23 Logsdon, 2002) and suggest that other factors, mainly changes in pore size distribution and 24 structure, may play a substantial role in specific situations. Indeed, ambivalent phenomena can 25 intervene simultaneously, which makes the understanding of the effective hydraulic properties 26 of stony soils difficult. The reduced volume available for flow might be partially compensated 27 by others factors. One compensation factor might be, as pointed out by Ravina and Magier 28 (1984), the creation of large pores in the rock fragments' vicinity. Indeed, the creation of new 29 voids at the stone-fine earth interface could generate preferential flows and hence increase the 30 saturated hydraulic conductivity (Zhou et al., 2009;Cousin et al., 2003;Ravina and Magier, 31 1984;Sauer and Logsdon, 2002). 32 hydraulic resistance of the stony fraction considering shape, size and orientation of inclusions 1 (the recommended value is 1.32 for clay soils according to ). 2 Two major characteristics are widely used to describe the hydraulic properties of unsaturated 3 soil: the water retention curve (ℎ) and the hydraulic conductivity curve (ℎ). These are both 4 non-linear functions of the pressure head h. One of the most commonly used analytical models 5 has been introduced by van Genuchten (1980), based on the pore-bundle model of Mualem 6 (1976), and given by: 7 In which ℎ is the pressure head were used: rock fragments (granite) with a diameter between 1 and 2 cm (1) and spherical glass 24 beads with a diameter of 1 cm (2) (see fig 1). The fine earth is classified as a clay (sand: 26%, 25 silt: 19%, clay: 55%). 26 Before each measurement campaign, fine earth was first oven dried for 24 hours at 105°C and 27 passed through a 2-mm sieve. To prepare a sample without any inclusion, fine earth was 28 compacted layer-by-layer to get an overall bulk density of 1.51 g/cm 3 (equal to the mean bulk 1 density of the fine earth measured in situ (Pichault, 2015)). For samples containing rock 2 fragments, stones were divided over four layers of soil application and laid on the fine earth 3 bed on their flattest side. The samples were then compacted layer-by-layer in a way that 4 maintains the same bulk density of fine earth as for samples without inclusions (as a result, the 5 global bulk density of samples varies according to stoniness). Even though the filling and 6 compaction procedure was conducted with precision, it is probably impossible to avoid local 7 bulk density heterogeneity as stones can move and/or soil between stones can be less 8 compacted due to difficult access of the area close to the stone during compaction. The same 9 procedure was to prepare samples containing glass balls. Once the specimen was made, it was 10 placed in a basket containing a thin layer of water during at least 24 hours in order to saturate 11 the soil from below. 12

Setup Description 14
We used the evaporation method to determine the unsaturated hydraulic conductivity and the 15 retention curve of a soil sample. The principle of this method is to simultaneously measure the 16 matric head at different depths and the water content of an initially saturated soil sample 17 submitted to evaporation. The total water loss as a function of time was monitored by a balance (OHAUS) with a 1 sensitivity of 0.2 g with an accuracy of ± 1 g with a time resolution of 15 min. A 50 W infrared 2 lamp was positioned 1 m above the sample surface to slightly speed up the evaporation process. 3 The light was turned off for the first 24 hours of every experiment, as the evaporation rate is 4 already high in a saturated sample. A measuring campaign lasted until three of the four 5 tensiometers ran dry (the tension sharply drops down to approximately a null value). At the end 6 of the experiment, the sample was oven dried for 24 hours at 105°C to estimate the . 7

Data Processing 8
A simplified Wind's method (1968) was used to transform matric potential and total weight 9 data over time into the hydraulic conductivity curve (Schindler, 1980 cited by Schindler and 10 Müller, 2006;Schindler et al., 2010). The method is further adapted in order to take into 11 account the data from four tensiometers. The method assumes that the distribution of water 12 tension and water content is linear through the soil column. It further linearizes the water 13 tension and the mass changes over time. The time step chosen to process the data is one hour. 14 By calculating the hydraulic conductivity based on measurements of two tensiometers and 15 linking it to the corresponding mean matric head, one can evaluate a point of the hydraulic 16 conductivity curve. We used every possible combination of two tensiometers (six here) to 17 obtain data points for the hydraulic conductivity curve. 18 Points of the hydraulic conductivity curve obtained at very small hydraulic gradients (defined 19 here as ∇ = ∆|ℎ| ∆ − 1) were rejected, because large errors occur in the near-saturation zone due 20 to uncertainties in estimating small hydraulic gradients (Peters and Durner, 2008;Wendroth, 21 1993). This highlights in its turn the necessity of reliable tensiometers to estimate the near- Wendroth, 1993). Using the least restrictive filter criterion (hydraulic gradient > 0.2) requires 25 fine calibration and outstanding performance of the tensiometers. Choosing a more restrictive 26 criterion leads to a larger loss of conductivity points, but provides more reliable and robust 27 data. We decided to use a filter criterion that does not consider hydraulic conductivity points 28 higher than the evaporation rate (from 0.1 to 0.2 cm/day in this case), resulting in a lower limit 29 of 1 cm/cm for the hydraulic gradient. 30 As pointed out by Wendroth (1993) and Peters and Durner (2008), the main drawback 31 associated with the evaporation experiment is that no estimates of conductivity in the wet range 32 can be obtained due to the typically small hydraulic gradients so that additional measurements 1 of the should be provided. To do so, we used constant-head permeability experiments (see 2 below). Except for the which is fixed using results from the constant-head permeability 3 experiments, the parameters of the VGM-model (1980) (Eq. (7)) are obtained by fitting 4 evaluation points from each combination of tensiometers using the so-called "integral method" 5 (Peters and Durner, 2006). 6 The soil sample used for permeability tests has the same size as the one from the evaporation 14 experiment (height: 65mm, diameter: 142 mm). A 2 cm thick layer of water was maintained on 15 top of the sample thanks to a Mariotte bottle. Water was collected through a funnel in a burette 16 and the volume of discharge was deduced from measurements after 30 and 210 min after the 17 beginning of the experiment (Δ = 180 min). 18

Numerical simulations 19
The HYDRUS-2D software was used to simulate water flow in variably saturated porous stony 20 soils. HYDRUS 2D solves the two dimensional Richards equation using the Galerkin finite 21 element method. 22 All the performed simulations assumed that rock fragments were non-porous so that "no-flux" 23 boundaries conditions were specified along the stones limits. Since we mimic the laboratory 24 setup, rock fragments were modelled as circular inclusions. The soil domain over which 25 simulations were performed had the same dimensions as the longitudinal section of the 26 sampling ring used in the laboratory experiments (14 x 6.5 cm). We considered the 2D fraction 27 of stoniness equal to the volumetric fraction. The parameters of fine earth used in the 28 simulations come from the fitting of the hydraulic conductivity and water retention curves 1 obtained in our laboratory experiments on stone-free samples (Table 1).  2 As a general rule, the hydraulic conductivity of a heterogeneous medium tends to be higher for 3 3D than for 2D simulations (Dagan, 1993). Similarly, for a same level of heterogeneity, the 4 flow will be more hampered using 1D rather than 2D simulations. In the present study, we 5 performed 2D simulations: the quantitative and qualitative conclusions from these experiments 6 can be only extended to the third dimension for their corresponding 3D form with an infinitely 7 long axis. However, we did keep the observation nodes on the edges and the larger sample size for the 19 following reasons. Firstly, we observed more changes in hydraulic gradient near stones. As 20 small variations of the hydraulic gradient can lead to substantial changes in the hydraulic 21 conductivity estimates, we chose to place observation nodes out of the influence of one specific 22 inclusion. This difficulty, especially at high stone contents, is the reason why the nodes are not 23 situated inside of the sample volume, but at the edges. Secondly, we checked whether the 24 pressure head was linearly distributed across the soil profile, which was the case. Finally, as we 25 are studying clayey soils, and as we are considering a pressure head range between pF 1.5 ad 26 2.5, these assumptions are likely to be fair enough (Peters et al., 2015). 27 As the relative mass balance error was large at the beginning of the simulations, we only started 28 considering values from the moment when this relative error was lower than 5%. This 29 validation criterion was set arbitrarily, based on the comparison between evaluation points from 30 the simulation of the evaporation experiment on stone-free samples and the expected values 31 obtained from the inputs of the simulation. The hydraulic conductivity curve was obtained 32 fitting the discrete conductivity data plus the simulated saturated hydraulic conductivity using 1 the so-called "integral method" (Peters and Durner, 2006), just like we did for the laboratory 2 experiment. 3

4
The was determined using a numerical constant-head permeability simulation. We 5 simulated a steady-state water flow of a saturated soil profile, with a constant head of 10 cm 6 applied on the upper boundary. The bottom boundary of the column was defined as a "seepage 7 face", which means that water starts flowing out as soon as the soil at the boundary reaches 8 saturation. The calculation method applied to the output data was identical to the permeability 9 experiment. 10 Table 2  and practical constraints in the numerical simulation, we added an increasing to observe the 20 evolution of the hydraulic conductivity curve. Simulations were performed on soil samples 21 containing 12 regularly distributed stones. One can notice that no investigations of the 22 unsaturated properties with coarse fragments above 30% of were performed. Indeed, given 23 that small variations of the hydraulic gradient can lead to substantial changes in the hydraulic 24 conductivity estimates, the tensiometers should be ideally positioned out of the direct influence 25 of one particular stone in order to obtain generalizable results. This implies the need for 26 relatively low stone contents (< 30% according to Zimmerman and Bodvarsson (1995)). 27

Treatments 11
Then, to study the relationship between saturated hydraulic conductivity, , and , we 28 performed five replications of four volumetric stone fractions (0, 20, 40 and 60%) with rock 29 fragments (1). We also tested a second type of inclusions, glass spheres (2), with a , of 20% 30 (1 replication). The first setup with rock fragments was concomitant with the one with glass spheres. Then, the four supplementary replications with rock fragments were processed for the 1 different volumetric fractions altogether: between replications the soil was oven dried for 24 2 hours at 105°C and passed through a 2-mm sieve. Numerical permeability simulations were 3 also performed involving 12 circular regularly distributed inclusions for the same (0, 20, 40, 4 60%). 5 Finally, we used supplementary numerical simulations to investigate the effect of the inclusion 6 shape and size on . To do so, simulations of the permeability test were performed on soil 7 containing stones of five different shapes: circular, upward equilateral triangle, downward 8 equilateral triangle, rectangle on its shortest side (L x 1.5L) and rectangle on its longest side 9 (1.5L x L)) with an of 10, 20 and 30%. We first performed simulations on soil containing 10 only one centered inclusion. We also performed permeability simulations on soil containing 12 11 and 27 regularly distributed inclusions (for each ). 12

Results and Discussion 13
In the following, results from laboratory experiments and numerical simulations will be 14 compared to the predictions of the different models presented in Section 2.1. The will be 15 represented by the median value predicted by the five models linking the properties of fine 16 earth to the ones of stony soil (Eq. (1) to Eq. (5)). This will be referred to as "results from the 17 predictive models" in the following and will be graphically represented by dotted lines. 18 The same predictive models assume that the shape parameters of the VGM-equations, n, l and 19 , do not depend on the stoniness, as suggested by Hlaváčiková and Novák (2014). As 20 mentioned above, unsaturated functions of stony soils have been barely studied. We will 21 compare results from unsaturated experiments and numerical simulations to predictive models 22 results following this assumption. 23 by Ravina and Magier (1984). Besides, we have to keep in mind that these elements are very 23 likely to have a different impact depending on soil texture, which was clay for both studies. 24

Effect of Stones on Saturated Hydraulic Conductivity 24
Glass beads were used to check the influence of rock characteristics on our conclusions about 25 . Since results with glass beads show a trend similar to the five replications with rock 26 fragments, we infer that it is not the rock fragment itself that produces bigger , but the 27 presence of a certain volume of inclusions. Besides, the variation observed between the trends 28 of the curves with rock fragments and glass beads could be due to the inner variation of the 29 hydraulic properties of samples, but it could suggest as well that depends on the shape and 30 the roughness of the inclusions. Nevertheless, we can only see the combined effect of these 31 factors in this experiment. This leaves the understanding of the major drivers of the and 32 their relative importance unclear. These elements are further investigated through numerical 1

simulations. 2
Besides the observed increase of with rock content, we can also observe a decrease in 3 between replications (see fig 3). In fact, as mentioned above, the global trend of increasing 4 is observed for each replication individually, but sampling procedure seems to have a large 5 impact on results too. There are significant differences (p<0.05) between replication 2 and 6 replication 5, the last one presenting lower . The drying and wetting cycles and/or the 7 sieving influence the hydrodynamic behavior of soil fraction since the effect decreases when 8 increases. This underlines the effect of soil texture and is an important aspect to take into 9 account in future studies. 10

Effect of the Stone Size and Shape on the Saturated Hydraulic Conductivity 11
To investigate the effect of the size of the inclusions and their shape on separately from 12 other factors of variation, we performed constant-head permeability simulations on samples 13 containing 1, 12 and 27 inclusions of various shapes, for a of 10, 20 and 30%. Table 3  14 illustrates the tendency of the effects and their respective factors. 15 Table 3  the higher the resistance to flow at a given stoniness. We suggest the decrease of is due to a 29 combination of the two following phenomena. The first one is the overlapping of the influence 30 zone of each inclusion, causing further reduction of . The concept of overlapping influence 31 zones was first proposed by Peck and Watson (1979) to explain higher decrease of the 32 hydraulic conductivity of stones very close to each other in comparison to isotropically 1 distributed stones. The second phenomenon could be that, for a given , the contact area 2 between stones and fine earth is higher for small stones than for bigger ones. Hence, a higher 3 tortuosity can be responsible for a lower flow rate. 4 The shape of the inclusions also has a visible impact on . For a fixed number of inclusions, 5 the is higher with rectangular inclusions on their shortest side and smaller with rectangular 6 inclusions on their longest side. Circular inclusions provoke a smaller reduction than triangular 7 inclusions. The orientation of the triangles does not have a pronounced effect on . Here 8 again, we observe a stronger effect of the size for higher stoniness. As an illustration, the 9 decrease in between circular and triangular inclusions is limited to 5% for a of 10% but 10 rises up to 14% for a of 30%. A similar behavior is observed with simulations including 11 either 1 or 27 fragments. 12 Considering a fixed of 20% (see Table 3 for samples containing only one spherical inclusion ( Table 3). The predicted 23 by the models is always higher than the determined by the simulations, except for soils 24 containing one inclusion on its shortest side. This can be a side effect of 2D simulations versus 25 3D measurements. Nevertheless, the numerical simulations show that the shape and the size of 26 inclusions may have an effect on , which is usually neglected by the current predictive 27 models. In general there is a concordance between models and simulations, whatever shape and 28 orientation of stones. This strengthens our hypothesis that macropore creation or heterogeneity 29 of bulk density close to the stones can occur and influence K se . Indeed, numerical simulations 30 cannot simulate the creation of voids, unless we create them manually and subjectively in the 31 domain. 32 Eventually, we hypothesize that, from a certain R v onwardsthe exact R v value depending on 1 the sampling procedure, the shape and roughness of inclusions, as well as soil texture -2 stoniness is at the origin of a modification of pore size distributions and of a more continuous 3 macropore system at the stone interface. This macropore system could overcome the other 4 drivers reducing . 5 predicted by the models for the corresponding . The hydraulic conductivity curves from the 9 predictive models and from the numerical simulations match hydraulic conductivity decreases 10 for increasing . According to these simulations, hydraulic conductivity in the unsaturated 11 zone is well defined using a correct and shape parameters do not depend on the stoniness. 12

Effect of Stones on Unsaturated Hydraulic
But this is not surprising since predictive models and numerical simulations rely on same 13 assumptions, i.e imperviousness of stones and an identical porosity distribution of fine earth. 14 As a result, these elements do not prove that shape parameters do not depend on the stoniness. 15 flattened than the ones measured on stone-free samples. This suggests that stones decrease 20 unsaturated hydraulic conductivity. However, it must be noted that we do not have unsaturated 21 K data for higher stone contents, whereas for , the effect of stoniness becomes more 22 obvious for > 20%. It might therefore be needed to find a way to conduct evaporation 23 experiments for higher stone contents in order to draw final conclusions. 24 In the numerical simulations, the presence of stones reduces the hydraulic conductivity in the 25 same way as predicted by the models, whatever the suction was. Similarly, the laboratory 26 experiments suggest that stones reduce the unsaturated hydraulic conductivity while laboratory 27 experiments in saturated conditions indicated that stones content might increase the . These 28 elements support the hypothesis of the macropore creation: according to the well-known law of could hypothesize that stones are always expected to decrease the hydraulic conductivity at low 1 effective saturation states. However, under saturated conditions, relationship between and 2 seems to be less trivial and requires further investigations considering soil texture and stone 3 characteristics. 4 5

Conclusion 6
Determining the effect of rock fragments on soil hydraulic properties is a major issue in soil 7 physics and in the study of fluxes in soil-plant-atmosphere systems in general. Several models 8 aim at linking the hydraulic properties of fine earth to those of stony soil. Many of them assume 9 that the only effect of stones is to reduce the volume available for water flow. We tested the 10 validity of such models with various complementary experiments. 11 Our results suggest that considering that stones only reduce the volume available for water flow 12 may be ill-founded. First, we observed that, contradictory to the predictive models, the 13 saturated hydraulic conductivity of the clayey soil of this study increases with stone content. 14 Besides, we pointed out several other potential drivers influencing , which are not 15 considered by these predictive models. We observed that, for a given stoniness, the 16 resistance to flow is higher for smaller inclusions than for bigger ones. We explain this 17 tendency by an overlapping of the influence zones of each stone combined with a higher 18 tortuosity of the flow path. We also pointed out the shape of stones as a factor affecting the 19 hydraulic conductivity of the soil. We showed that the effect of the shape depends on the 20 inclusion size and inversely that the effect of inclusion size depends on its shape. Finally, our 21 results converge to the assumption that this contradictory variation of could find its origin 22 at the creation of voids at the stone-fine earth interface as pointed out by Ravina and Magier 23 (1984). Even if the very mechanisms behind these observations remains unclear, they seem to 24 strongly depend on , shape and roughness of inclusions. However, as we conducted these 25 experiments on a specific clay soil only, and given the fact that structural modifications are 26 textural dependent, our results can't be extrapolated to other soil textures without similar 27 experiments. Finally, as we worked with disturbed samples, our results do not include 28 quantification of natural phenomenon such as swelling and shrinking that occurs naturally for 29 clay soils. 30 These findings suggest that the aforementioned predictive models are not appropriate in all 31 cases, particularly under saturated conditions. Models should take into account the required in order to explore the hydraulic properties of stony soils and to develop new models 1 or adapt the existing ones. The direct observation of undisturbed stony samples porosity using 2 X-ray computed tomography or magnetic resonance imaging could confirm -at first -and then 3 help better understand the mechanism of supposed voids creation at the stone-fine earth 4 interface. However, under unsaturated conditions, these considerations should be more