Root induced changes of effective 1D hydraulic properties in a soil column

Roots are essential drivers of soil structure and pore formation. This study aimed at quantifying root induced changes of the pore size distribution (PSD). The focus was on the extent of clogging vs. formation of pores during active root growth. Parameters of Kosugi’s lognormal PSD model were determined by inverse estimation in a column experiment with two cover crops (mustard, rye) and an unplanted control. Pore dynamics were described using a convection–dispersion like pore evolution model. Rooted treatments showed a wider range of pore radii with increasing volumes of large macropores >500 μm and micropores <2.5 μm, while fine macropores, mesopores and larger micropores decreased. The non-rooted control showed narrowing of the PSD and reduced porosity over all radius classes. The pore evolution model accurately described root induced changes, while structure degradation in the non-rooted control was not captured properly. Our study demonstrated significant short term root effects with heterogenization of the pore system as dominant process of root induced structure formation. Pore clogging is suggested as a partial cause for reduced pore volume. The important change in micro- and large macropores however indicates that multiple mechanic and biochemical processes are involved in root-pore interactions.

Abbreviations PSD pore size distribution h m median pressure head r m median pore radius σ standard deviation of PSD θ s saturation water content θ r residual water content K s saturated hydraulic conductivity V drift term λ dispersivity

Introduction
Soil hydraulic properties are the common result of particle size distribution (texture) and aggregation (structure). Beside texture, soil structure is the main property to shape the fundamental relations for soil water flow, i.e. retention and hydraulic conductivity, in the saturated and near-saturated range (Cresswell et al. 1992). Among the various driving factors of soil structural porosity, vegetation is playing a dominant role. Roots are a key element in plant related effects on soil structure and soil hydrology (Gregory 2006;Bengough 2012;Logsdon 2013). Recently, Carminati et al. (2010) and Moradi et al. (2011) demonstrated that contrary to the usual assumption of higher depletion in vicinity of roots (Gardner 1960), soil water content is higher in rhizosphere soil compared to bulk soil over a wide range of pressure heads, indicating a significant change of the water retention curve in vicinity of plant roots. Several pathways of root influence on soil hydraulic properties have been described. Temporal pore clogging occurs due to roots growing into pre-existing pores (e.g. Gish and Jury 1983;Morgan et al. 1995). Scanlan (2009) suggested that root in-growth results in the division of larger into smaller pores. Micro-fissures and cracks are structural pores formed in the root zone of transpiring plants by more intense wetting-drying (e.g. Dexter 1987;Mitchell et al. 1995;Young 1998;Whalley et al. 2005). These effects are depending on the lifespan of roots. Meek et al. (1990) and Murphy et al. (1993) measured reduced infiltration rates as long as plants are actively growing and their roots block pore channels. After root decay, bio-macropores and root-induced micropores are formed (Cresswell and Kirkegaard 1995;Mitchell et al. 1995;Wuest 2001;Horn and Smucker 2005;Ghestem et al. 2011). These pores have high connectivity (Pagliai and De Nobili 1993;Whalley et al. 2005), thereby facilitating water transport through the soil (Gish and Jury 1983;Murphy et al. 1993;Suwardji and Eberbach 1998). As a result, e.g. Disparte (1987) measured an increase of infiltration with higher rooting density after decomposition of roots. Gyssels et al. (2005) confirmed this finding by establishing a direct relation between root density and reduction of soil erosion.
The extent of root induced changes of soil hydraulic properties is influenced by both soil and root characteristics. Scanlan (2009) did not find any root effect on soil hydraulic properties in a column experiment using a sandy substrate. We suppose that changes of pore properties depend strongly on (i) the extent of existing growth paths for root penetration (Feeney et al. 2006) and (ii) the relation between root volume and pore volume (Bengough 2012). Yunusa and Newton (2003) reported differences among species in their effects on soil hydraulic properties. Perennials and woody plants substantially changed flow behavior while annual crops had hardly any influence. Among annual plants they suggested root diameter as main trait for effectively priming soil hydraulic properties. Higher strength of coarse roots allows more effective shift of soil particles and lower tendency of root buckling under mechanical stress (Clark et al. 2003). Coarser root axes exert higher radial pressures among soil penetration, thereby enlarging existing pores while increasing density of adjacent rhizosphere soil (Dexter 1987;Archer et al. 2002;Kirby and Bengough 2002;Whalley et al. 2004).
A major challenge for understanding root effects on soil hydraulic properties is measurement. Under field conditions other environmental and management effects can mask the distinct influence exerted by plant roots (Bodner et al. 2013). Therefore disturbed soil is often preferred to undisturbed field samples to observe root induced changes under controlled conditions (e.g. Scanlan 2009). In spite of more powerful 3D imaging methods to capture small scale pore processes in recent years (e.g. Gregory et al. 2003;Mooney et al. 2012), still scaling to macroscopic parameters governing the effective hydraulic behavior of soil is challenging. At the macroscopic scale inverse methods have been increasingly used to characterization hydraulic properties from observed state variables in a soil sample (e.g. Hopmans and Šimůnek 1999;Hopmans et al. 2002;Ritter et al. 2003). In spite of being indirect estimates only, these effective parameters can properly describe the hydraulic material properties (water retention, hydraulic conductivity) that underlie an observed the flow behavior. Comparison to direct measurements (pressure plate) showed that they deviated from inverse retention curves mainly towards the dry end (Šimůnek et al. 1998). In the wet range, where soil structural dynamics have key role, there was less deviation. Reliability of inverse estimates also depends on the range of data used for optimization estimation. For example, tension infiltrometers cover a narrow near-saturated pressure head range only (h>−20 to −15 cm; e.g. Schwen et al. 2011). Drainage experiments from large columns typically achieve an intermediate range (h>−150 to −100 cm; e.g. Ritter et al. 2004).
Outflow experiments from small sample cylinders can cover a relatively wide pressure head range (h > −800 cm; e.g. Eching and Hopmans 1993). Šimůnek et al. (1998) showed that inverse parameters could better describe flow processes in the natural soil systems compared to lab samples. Furthermore inverse parameters are often obtained from sampling approaches that cover a comparatively high representative elementary volume. This is an important advantage for representative characterization of the highly variable structural range.
From an agricultural point of view biological management of soil structure by roots ("biotilling") is still at its infancy. Yunusa and Newton (2003) presented the concept of primer-plants, i.e. plants without a direct economic benefit, but effective in conditioning the soil for cash crops and in conserving environmental resources. Currently this type of plants is used as cover crops in agro-environmental programs to minimize nitrate leaching (e.g. Vidal and López 2005) and reduce soil erosion (e.g. Zuazo and Pleguezuelo 2008). Several authors observed cover crop effects on soil structural properties such as aggregate size and stability (Liu et al. 2005) as well as hydraulic processes such as water infiltration (Carof et al. 2007;Bodner et al. 2008). Williams and Weil (2004) showed that cover crops could be an effective way to alleviate soil compaction due to root biopores being used by the following soybean crop to penetrate the soil.
The objective of our study was (i) to identify if there is a short term influence of growing plants on soil hydraulic properties in a soil column experiment, (ii) to analyse which particular pore property and pore range are influenced during the phase of active root growth, and (iii) to model the temporal change of the pore size distribution (PSD) in a non-planted compared to a rooted soil column. Our model plants were two commonly used cover crop species. The focus of the study was on the extent of pore clogging vs. formation of new pores by actively growing roots. We expected higher changes in soil hydraulic properties with more densely rooted soil, and a dominant reduction of pore volume in the pore range corresponding to root diameters. The overall aim was to improve the quantitative understanding of hydrological interactions of roots with soil structural porosity.

Column experiment
The influence of plant roots on soil pore properties was measured in a column experiment with planted vs. unplanted soil columns using an inverse evaluation procedure to quantify hydraulic properties.

Soil column design
We used twelve cylindrical (60 cm high, 15 cm diameter), custom-built soil columns made of Plexiglas. The column design (Fig. 1) was similar to the one used by Kosugi and Inoue (2002), Ritter et al. (2004) and Yang et al. (2004).
The columns were equipped with TDR probes and mini-tensiometers (LOM, EasyTest, Poland) in six Each column (Ø 15 cm, height 60 cm) is equipped with six TDR sensors and six tensiometers (LOM, Easy Test) in 10 cm increments, a suction plate (−800 cm air entry) at bottom and microdrip-irrigation at the top. (Manufacturer of columns, suction plate and drip irrigation: Technisches Büro für Bodenkultur, Austria) depths separated by 10 cm increments (5,15,25,35,45 and 55 cm) for continuous measurement of water content and matrix potential. At the top of each column a microdrip-irrigation was installed to apply a controlled upper boundary flux for drainage experiments and to supply plants with water during the growing period. Irrigation was managed by a PC-controlled multi-channel precision pump with each column supplied by a separate channel. Growth lamps (4 × Narva LT 58 Watt,8 × Sylvania GroLux F58 Watt) were mounted above the columns to lighten plants. A simple timer controlled lightening duration which was set at 10 h per day.
For controlled drainage, each column was equipped with a porous ceramic suction plate at the bottom to apply a defined lower boundary pressure head condition via different setting of the connected vacuum pump (0 to −800 mbar). From the suction plate a pipe led to a collection bottle (volume: 2,000 cm 3 ) placed on a balance for recording the outflow volume in 10 min intervals. Pipes were fixed to a frame in order to avoid unmeant pressure on the balance plates that could bias outflow measurement.
All components of the system (TDR, tensiometers, suction at bottom boundary) were previously calibrated and tested. Products and systems of the various manufacturers (in total 144 sensors) were combined via a selfprogrammed software based on MS.NET 2.0 for data collection at 10 min intervals. The dataset was finally stored in an MS Access database.

Experimental conditions
All 12 soil columns were filled with air dried soil sieved to <2 mm. The soil type was a calcareous chernozem with a silt loam texture (0.19 kg kg −1 sand, 0.56 kg kg −1 silt, 0.24 kg kg −1 clay). Organic carbon content of the soil was 0.025 kg kg −1 . In order to exclude complex interactions between roots and soil micro-organisms, the experiment was performed with sterilized soil. Heat sterilization was done using an autoclave. Preliminary tests did not show differences in early vigour of plants growing in autoclaved and non autoclaved soil. Soil sterilization ensured that root effects were limited to mechanical influences (pore clogging, pore enlargement by root pressure, entanglement of soil particles) and root-induced gluing of soil particles by root exudates.
Column filling was performed carefully to avoid displacement of fine soil particles and layering. A constant bulk density of 1.3 g cm −3 was adjusted by successive filling and compaction of 5 cm layers. Boundary compaction was minimized by loosening the contact zone at each filling level. Before sensor installation, columns were fully saturated and drained for three times to limit further soil settlement during the main experimental run. The top of the columns was covered with a two cm layer of fine gravel to avoid surface effects from drip irrigation as well as to reduce evaporation.
After an initial drainage experiment with unplanted columns, mustard (Sinapis alba L.) and rye (Secale cereale L.) were planted in four replicates, while four columns remained unplanted. These two commonly used cover crops differ in root architecture, mustard having a taproot system and rye a fibrous root system. Plant density followed common seeding rates as used also in a parallel field experiment (Bodner et al. 2013). Three vigorous plants per column were maintained after an initial seeding of six seeds per column. After 3.5 months roots could be observed at the bottom of the columns indicating a sufficiently developed root system. At this time, rye was at BBCH stage 32-33 and mustard at the onset of flowering (BBCH 60-61). After harvest of shoot biomass, a final drainage experiment was performed.

Measurements
Soil hydraulic properties were determined in a drainage experiment following Kosugi and Inoue (2002) and Ritter et al. (2004). An initial measurement (without root effects) was performed before planting, while the final measurement (with root induced changes) was done after full crop development. The lower boundary during the drainage experiment was set to a constant pressure head of −500 cm. At the upper boundary an irrigation flux of 5 mm h −1 was applied until a steady outflow rate was achieved as initial condition. Thereafter irrigation was stopped and water content, pressure head and outflow rate were monitored during the redistribution phase. Soil surface was covered by parafilm to avoid evaporation (i.e. no flux upper boundary condition), while at the lower boundary drainage continued under a constant pressure head of −500 cm. The transient profiles of water content, pressure head and outflow during the redistribution phase entered the inverse evaluation procedure for hydraulic property estimation. Data recording was stopped after about 100-110 h when outflow rates had approached zero and changes in water content over a period of 6 h were as low as 0.001 cm 3 cm −3 . A second drainage experiment was performed subsequently to obtain a validation data set. Columns were again irrigated to the initial steady outflow rate and drained under the same conditions as described previously.
After the end of the experiment, soil was removed from the columns and root morphological parameters were measured at 10 cm increments. Root analysis followed the procedure described by Himmelbauer et al. (2004). Roots were washed free from soil over a set of sieves up to 0.5 mm with tap water, stained and thereafter analyzed using WinRhizo (Regent Instruments Inc.). Finally root and aboveground dry matter were measured after ovendrying at 60°C until constant weight.

Determination of hydraulic properties
Determination of soil hydraulic properties was done by inverse parameter estimation using HYDURS 1D v. 4.16 (Šimůnek et al. 2013). The procedure is described in detail by Hopmans et al. (2002). In short, parameters for the constitutive relations -soil water retention S e (h) and hydraulic conductivity K(h)in a 1D vertical water flow simulation via Richards' equations are obtained by minimizing the deviation between observed and simulated state variables. We used the model of Kosugi (1996) to describe S e (h) and K(h) which is based on a lognormal pore-size distribution. S e (h) is given by where S e (−) is the effective saturation corresponding to θ−θ r θ s −θ r with θ r (cm 3 cm −3 ) being residual water content and θ s (cm 3 cm −3 ) saturation water content. Erfc is the complementary error function, h m (cm) the median pressure head and σ (−) the standard deviation of the logtransformed pressure head. K(h) can be written as where K s (cm s −1 ) is saturated hydraulic conductivity and l (−) is a tortuosity factor.
Observed state variables for inverse optimization were the transient time series of water content, pressure head and cumulative outflow during the redistribution phase in six soil layers and 10 min interval. Simulated values of these variables were obtained by numerically solving Richards' equation in HYDURS 1D with a no flux upper boundary and a constant head (−500 cm) lower boundary condition. The total of differences between measured and simulated values are expressed by the objective function being where the right side represents residuals between measured (y j *) and corresponding predicted (y j ) variables using the soil hydraulic parameters of the optimized parameter vector β. Measured state variables are denoted by m y , whereas the number of measurements for a certain state variable are given by the variable n j . Weighing factors v j and w ij can be included for a given state variable or an individual data point respectively. For our study, the total number of data points in the objective function was between 7,600 and 7,900. The inverse problem was solved by minimizing the objective function using the nonlinear Levenberg-Marquardt optimization algorithm implemented into HYDRUS 1D. HYDRUS 1D allows a maximum number of 15 parameters to be optimized simultaneously. Still it is recommended to reduce the number of parameters to avoid non-unique solutions (e.g. Durner et al. 1999;Abbasi et al. 2003a, b). Therefore we subdivided the soil column into two layers only (0-30 cm, 30-60 cm), assuming distinct temporal changes between the more densely rooted upper part and the less rooted lower part of the columns. Residual water content was set to 0.067 and the tortuosity parameter l was set to 0.5 for all simulations as predicted by the texture based pedotransfer function Rosetta (Schaap et al. 2001). In this way a total number of eight parameters (θ s , h m , σ and K s for two layers) had to be optimized. Initial parameter estimates were obtained from direct fitting of the Kosugi retention model to data-pairs of water content and pressure head using RETC (van Genuchten et al. 1991). Initial estimates of K s were also taken from Rosetta.
In order to ensure a unique solution (global minimum) we changed the initial values and assessed whether the final estimates converged at the same values.
Furthermore we checked the statistical parameters (standard error coefficient, confidence limits, R 2 between measured vs. predicted) provided by HYDRUS 1D to ensure reliability of estimation results. Finally also statistical analysis of the obtained parameter values provides an indicator for the quality of optimization results. In case of a non-unique solution, error variance would increase by an additional random effect due to different optimization quality, while the relative weight of the fixed factor effects (e.g. plant influence) would decrease. Thus also statistical significance indirectly reveals the reliability of the parameter optimization results.

Simulation of pore evolution
The pore evolution model presented by Or et al. (2000) was tested to describe the root induced changes of hydraulic properties. Pore size distribution (PSD) f is the first derivative of the retention curve and can be written as where r (μm) is the pore radius, r m (μm) is the median pore radius, and σ (−) is its standard deviation. The median pore radius (r m ) can be calculated from the median pressure head h m using the Young-Laplace equation.
According to the pore evolution model of Or et al. (2000) the temporal change of the PSD is given by where t is time, V (μm s −1 ) is a drift term, D (μm 2 s −1 ) a dispersion term and M (s −1 ) a degradation term. The drift and dispersion terms quantify changes with time of r m and σ respectively. M represents a sink term for changes in total porosity. Dispersion is related to drift by a constant dispersivity λ (μm). An analytical solution of Eq. 5 was derived by Leij et al. (2002). To model the distinct temporal changes of the PSD between the initial and final state of the three treatments (no plant, mustard, rye) we optimized the cumulative drift term T (i.e. the integral of V; cf. Leij et al. 2002) as well as dispersivity. M was set equal the reduction in total porosity. While other authors limited degradation to the macropore range (e.g. Schwärzel et al. 2011), due to the lack of proper data, we did not attribute degradation to any distinct pore range. After testing the model for the temporal changes between the initial and final state, we simulated the case of pore evolution with increasing rooting density based on the data from the final experiment. For this purpose we used a logistic function following Leij et al. (2002) to describe the drift term as a function of root length density while dispersivity was again assumed as constant over time. All calculations of pore evolution were done with Matlab Version 8 R2012b.

Statistical evaluation
Statistical data evaluation (hydraulic parameters, plant data) was performed by analysis of variance with the procedure PROC MIXED in the software SAS 9.2. This procedure is based on restricted maximum likelihood estimates of the variance components (Littell et al. 1998) and provides Wald-type F-statistics using GLSE (generalized least squares). A mixed model for performing a proper analysis of variance is required as our data include repeated measures over depth (two layers) and time (initial and final experimental run). Thus an adequate correlation model has to be fit to account for serial correlation of non-randomized repeated measurements on the same experimental unit (Piepho et al. 2004). According to the Akaike Information Criterion (AIC) an unstructured (UN) model fitted best to our data. In case of significant effects at p<0.05 in the analysis of variance, comparison of means was performed using a two-sided t-test.

Plant development
Shoot growth of the plants grown in the columns achieved dry matter values of 9.19±1.03 g for mustard and 1.40±0.22 g for rye. This would correspond to 5,200 kg ha −1 dry matter for mustard and 792 kg ha −1 for rye, being in the range of field measured values (e.g. Bodner et al. 2010Bodner et al. , 2013. While mustard plants had a vigorous aboveground growth, tillering of rye was less compared to field grown plants, resulting in a lower biomass than expected at the end of the experiment (BBCH 32-33). Figure 2 shows root length density of the two species at the end of the experiment. Root length density of mustard was significantly higher compared to rye in the upper soil layer. Visual inspection during root analysis showed that the extension of the mustard tap root was around 10 cm, i.e. within the uppermost soil layer. Generally both plants had a lower rooting density compared to plants grown under field conditions (e.g. Bodner et al. 2010).
Rye had a significantly higher root diameter than mustard in all depth. Root diameters of both species (mustard 0.27 mm, rye 0.34 mm) were in the range of values measured under field conditions (Bodner 2007). Within each species there was no significant decrease of root diameter over column depth.
Total root volume was 6.36±1.40 cm 3 for mustard and 0.79±0.02 cm 3 for rye. Scanlan and Hinz (2010) used a lognormal function to study species differences in root volume distribution over different diameter classes. Our sample showed a median root diameter of the lognormal distribution at 0.21 mm for mustard and 0.16 mm for rye, and a standard deviation of the distribution of 0.94 for mustard and 0.72 for rye respectively. The higher median diameter of mustard shows that in spite of the slightly smaller average diameter, coarse root segments, particularly the tap root, substantially contributed to total root volume. The higher standard deviation of mustard indicates a more even contribution of coarse and root fine axes to total root volume. The narrower distribution for rye reveals that here root volume was built by morphologically less differentiated axes types.
Estimation of soil hydraulic properties Table 1 shows the direct parameter estimates obtained by fitting the Kosugi retention model to data pairs of water content and pressure head using RETC.
These direct estimates were then used as starting values for inverse optimization of the Kosugi parameters in each soil column. Table 2 gives the obtained parameter values. Additionally Fig. 3a and b show measured and simulated outflow for the calibration and validation data sets.
With only some exceptions a high overall R 2 ( Table 2) between observed and predicted data was obtained. This is shown for cumulative outflow in Fig. 3a, b. Data of simulated and measured water content and pressure head time series in the six observation depths are not shown. The overall optimization R 2 did not show significant differences between the initial and final experiment in spite of an average higher R 2 for the final run. Importantly there was no influence of plant treatments on the optimization R 2 . Also the standard errors of optimized parameters were low. On average the relative standard error (i.e. mean/standard error) for θ s was 0.53 %, for h m 3.0 % and for σ 1.6 %. Only K s had a higher standard error (22.1 %) indicating highest uncertainty for this parameter.
For outflow profiles further statistical indicators for the quality of estimated vs. measured data were determined (Table 3) using the IRENE software (Fila et al. 2003). They also suggested a generally good fit between measured and simulated data and thus reliable hydraulic parameter estimates obtained by optimization. A slightly lower quality of fit was obtained for the validation compared to the calibration data, although this difference was statistically significant only for the relative mean error (4.0 % vs. 6.1 %). Again no indicator showed any significant difference between plant treatments. Thus fitting criteria demonstrated that optimization induced no systematic bias for treatment comparison.
The reduction in average total porosity between initial and final experiments (mean of all treatments at initial experiment: 0.45 cm 3 cm −3 , mean of all treatments at final experiment: 0.41 cm 3 cm −3 ; cf. Table 2) indicates soil settlement over time. Also mean h m decreased in the course of the experiment from −162.7 to −306.8 cm indicating a shift to a smaller pore size, while σ increased from an average of 2.10 at the beginning to 2.65 at the end. Interestingly there was a significant increase in K s between the initial and the final experiments in spite of a reduced θ s and h m . Among Kosugi pore parameters, pore radius standard deviation showed the strongest relation to K s (r 2 =0.57 and 0.50 for the initial and final experiment, respectively). Thus in spite of the decrease in h m and the related shift of median pore radius to smaller pore classes, the broader range of different pore radius classes indicates the formation of additional macroporosity. These pores of large diameter can essentially contributed to K s due to high   The validation data set. The calibration data were used to inversely estimate soil hydraulic property parameters, while the simulated outflow for the validation data was obtained in a forward simulation using the same hydraulic parameters. All simulations were done with HYDRUS 1D transport capacity (Poisseuille's law) and lower tortuosity (Vervoort and Cattle 2003).
Plant root influence on pore size distribution Table 4 shows the results of analysis of variance highlighting which factors significantly influenced the inversely determined hydraulic parameters. The corresponding parameter values are given in Table 2. A significant influence of plant treatments was found for σ and K s . Plant treatments had a significant interaction with time, i.e. differences were only relevant for the final experimental run (after root growth) as expected. For both parameters, treatments with plants (σ=2.96; K s =26.4 cm min −1 ) differed from the non-planted control (σ=2.02; K s =1.6 cm min −1 ), while planted columns were similar among each other. The temporal change of these parameters between the initial and final experiment was not significant for the unplanted columns (σ initial =2.32 vs. σ final =2.02; K s,initial =1.1 cm min −1 vs. K s,final =1.6 cm min −1 ), while they changed significantly for the planted ones (σ initial =2.00 vs. σ final = 2.96; K s,initial =6.2 cm min −1 vs. K s,final =26.4 cm min −1 ).
For saturation water content plant treatments showed significant interaction with time and depth. Comparison of means revealed that a plant effect could be demonstrated only for the final experiment and the more densely rooted upper layer. Again the planted treatments (mustard 0.43 cm 3 cm −3 , rye 0.41 cm 3 cm −3 ) differed significantly from the non-planted soil (0.39 cm 3 cm −3 ) while being similar among each other. For θ s temporal dynamics showed a significant reduction between the initial (0.47 cm 3 cm −3 ) and final experiment (0.40 cm 3 cm −3 ) in the lower layer for all treatments, while in the upper layer only the bare soil treatment showed a significant pore loss. The planted columns obviously reduced the decrease in θ s via roots stabilizing the pore system against further settlement.
H m was not significantly changed by plant roots and only differed in average between the initial (−162.7 cm) and final experiment (−306.8 cm) as well as between the two layers (h m,upper =165.3 cm vs. h m,lower =304.2 cm).
Changes in overall PSD between the initial and final experiment for the two depths (i.e. interaction of T*D) and the different treatments (i.e. interaction of T*P) are shown in Figs. 4 and 5. Tables 5 and 6 give the Calculations are based on the agreement between measured and modelled cumulative outflow for the calibration and validation data sets (RMSE is root mean square error and CD is coefficient of determination with 1 showing optimum fit) respective volumetric changes for different pore classes using the classification of the Soil Science Society of America (SSSA 2013). This classification defines macropores having a radius >37.5 μm (37.5≤ r < 500 μm very fine macropores), mesopores between 37.5 μm and 15 μm, and micropores <15 μm (r< 2.5 μm ultramicropores and cryptopores).
There was a general decrease in pore volume over time in both layers (upper: −0.013 cm 3 cm −3 ; lower −0.063 cm 3 cm −3 ). As shown in Table 5, this was the result of a reduction in a pore range between 2.5 and 500 μm (micropores to very fine macropores). An increase was registered for micropores <2.5 μm and macropores ≥500 μm. The lower layer showed 79 % higher overall loss of pores while increase in pores <2.5 μm and ≥500 μm was similar to the upper layer. Total pore loss in the lower layer was 13.5 % of initial porosity, while in the upper layer it was only 3.1 %.
There was a clear differentiation in pore dynamics between the planted and unplanted columns (Fig. 5, Table 6).
The reduction in soil porosity for the non-planted columns covered the whole range of pores. Total loss of porosity was −0.058 cm 3 cm −3 increasing towards the range of macropores and finer micropores <2.5 μm. The rooted columns on the contrary, showed a reduction in the range of larger micropores to very fine macropores (2.5 to 500 μm), while increasing in the finer micopore classes <2.5 μm as well as in larger macropores ≥500 μm. Reduction of larger micropores and mesopores was substantially higher in the planted in Pore size distribution (log-log scale) of the non-planted and planted (mustard, rye) soil columns at the initial (before planting) and final (with plant influence) state. Temporal changes in volumetric frequency at different pore radius classes are highlighted (black colour: decreasing frequency, light grey colour: increasing frequency) relation to the unplanted columns, while reduction in very finer macropores was similar. Mustard reduced total porosity to a higher extent (−0.032 cm 3 cm −3 ) compared to rye (−0.024 cm 3 cm −3 ). The main differentiation between the planted treatments was in the volumetric increase in larger macropores >500 μm. Here mustard had a 75 % higher increase than rye.

Simulation of root induced pore evolution
Optimization of the governing parameters in the pore evolution model of Or et al. (2000) resulted in a cumulative drift of 15.35 μm and a dispersivity of 1.91 μm for the non-planted treatment. For the planted treatments the obtained parameters were 4.20 and 0.38 μm for mustard and 2.45 and 0.42 μm for rye, respectively. Measured and predicted evolution of PSD is shown in Fig. 6. Although the simulated PSDs suggest a satisfactory prediction by the model, for some cases the approach was not appropriate. Particularly the non-planted treatment showed a reduction in σ over time in most cases, i.e. a tendency to a narrower PSD. This could not be reproduced by a model based on a dispersion like process of pore evolution, resulting in a high residual error in optimization and therefore unreliable parameter estimates.
Our measurements showed pore evolution towards a broader PSD in planted columns being expressed by an increase in σ. This differentiation was significant between non-rooted (σ=2.02) and rooted soil. The trend to higher σ in more densely rooted soil of mustard (σ=3.13) compared to rye (σ=2.78) was not significant. This dispersion like dynamics allowed application of the pore evolution model. Figure 7a shows measured and modelled changes in PSD between a non-rooted (RLD 0 =0 cm cm −3 ) and a rooted soil with increasing root length density. For this example we assumed that the distinct root influence of the two species was only due to their different rooting density (RLD 1 =0.14 cm cm −3 , RLD 2 =1.93 cm cm −3 ).
For this simulation we did not optimize the cumulative drift T, but we used a drift term V as a logistic function of root length density (Fig. 7b) following Leij et al. (2002). The sharp increase of the drift function to a maximum reveals that there was an immediate reduction in median pore radius between the non-rooted   Figure 7a demonstrates that there is an increasing deviation of the modelled PSD with higher root length density. This is the result of assuming a constant λ and thereby a functionally similar change of σ with RLD as used for the drift V. However measurements demonstrated that there was a trend towards higher s with increasing RLD between the two species. This required a distinctive λ for each predictive step which resulted in an improved model fit. This indicates that the linkage between V and D is problematic in case of root driven pore evolution because of qualitatively and/or quantitatively distinctive changes in r m and σ.

Discussion
Roots are a key factor in soil aggregation and formation of structural porosity (Angers and Caron 1998). The present study analyzed root influences on PSD during active root growth. The main objective was the quantitative assessment of changes in PSD induced by a tap rooted dicot (mustard) and a fibrously rooted monocot species (rye) compared to a non-planted control.

Plant development
Root development in the columns achieved lower density compared to field grown plants. This might be related to the unstressed conditions in the columns with optimum water supply and regular addition of nutrient solution during the experiment. De Willigen et al. (2000) reported that a rooting density between 0.5 and 1.0 cm cm −3 might be sufficient for optimum plant supply under non-limiting conditions. The moderate root development of rye, with root length density below 1.0 cm cm −3 in all layers, however is also related to low tillering and thus little development of shoot-born roots. Beside maximum rooting density, the two species also differed in average root diameter. Root effects on soil porosity are strongly related to root diameter. While coarser roots have higher strength to shift soil particles due to lower tendency of root buckling under mechanical stress (Bengough et al. 2011), finer lateral roots can access smaller sized pores. Using a lognormal distribution model of root volume over different root diameter classes (Scanlan and Hinz 2010) we highlighted the different volume allocation between the tap rooted and fibrous rooted species: a higher standard deviation of the distribution for mustard pointed to an even allocation of root volume to coarse (tap root, first order laterals) and fine axes (higher order laterals). The median root radius of mustard was 56 % higher than its average radius, while for rye there was hardly a difference between both values. This reveals the important contribution of the tap root to overall root volume for mustard. Morphological differentiation between root axes of rye was less (narrow standard deviation), i.e. the bulk of root volume was formed by axes of similar diameter. Beside the generally more homogeneous diameter distribution in monocot root systems, this distribution pattern also expressed the dominance of primary and basal axes with similar morphology (Kutschera et al. 2009). Coarser shoot-born From root characterization we could expect stronger overall impact on soil porosity by mustard due to higher rooting density. The lower average diameter and larger standard deviation of root volume distribution of mustard (allocation to fine axes) pointed to higher potential of root in-growth and pore clogging. Macropore effects could be expected from coarser rye axes (higher average diameter) as well as the tap root of mustard.

Soil hydraulic property determination
Parameters of the soil PSD were obtained by inverse estimation from a drainage experiment. Interestingly only for the initial experiment (unplanted columns) there was a significant relation between the directly determined starting values (Table 1) and the final inverse estimates (Table 2), while for the final experiment (planted columns) parameter values did not correlated. Direct parameter estimation from a drainage experiment with large soil columns involves substantial uncertainty due to extrapolation from a rather narrow, near saturated range of measured data towards the entire retention curve (Ritter et al. 2004). In our case the pressure head range was between 0 cm and a minimum at the lowest tensiometer (55 cm depth) of −235.4 cm for the initial experiment and −149.9 cm for the final experiment. The narrower range of θ-h data pairs in the final experiment might have decreased the reliability of direct estimation. Also enhanced soil structuring over time might have been captured appropriately only by the inverse procedure, explaining the lack of correlation between the two procedures at the final experiment.
Most work on inverse estimation of soil hydraulic properties has been performed on small sample cylinders (e.g. Durner et al. 1999;Durner and Iden 2011). Experiments in small samples (e.g. one-step, multi-step outflow) cover a wider range of θ and h compared to larger volumes (e.g. soil columns, lysimeters). Also equilibration of the system to an applied boundary pressure head is achieved more readily. The pressure head range reported from studies using soil columns for inverse parameter estimation was from saturation to about −100 cm (e.g. Kosugi and Inoue 2002;Ritter et al. 2004;Yang et al. 2004;Köhne et al. 2005;Javaux and Vanclooster 2006). Using a suction plate with a lower boundary pressure head of −500 cm we obtained a minimum pressure head between −149.9 and −235.4 cm. Thus in our study the range of water content and pressure head data was relatively wide for a column experiment. Large soil column experiments also face higher vertical heterogeneity compared to small cylinders (e.g. Inoue et al. 2000;Abbasi et al. 2003b;Ritter et al. 2004;Schwärzel et al. 2006). This increases the number of parameters to be estimated (seven to twelve in the cited studies; in our case eight).
A main concern in inverse modelling is convergence of the optimization algorithm at a global minimum. Particularly when the number of optimized parameters is high, the risk of non-unique solutions increases. Some authors therefore used global optimizing approaches to better explore the parameter space and avoid convergence at a local minimum (e.g. Vrugt et al. 2008). Still when properly controlling fit statistics as well as varying starting values, also local search algorithms such as the Levenberg-Marquardt algorithm implemented in HYDURS 1D and used in this study can provide appropriate parameter estimates . We observed the highest standard error for K s , while it was substantially lower for the retention parameters. Other studies with similar experimental setup did not report this indicator of parameter quality for comparison. Still overall R 2 and SSQ values of the minimized objective function were in the same range as those reported by Abbasi et al. (2003b) who optimized up to 12 parameters. Other studies with different setup also reported highest uncertainty for K s among optimized parameters (e.g. Vrugt et al. 2003). We suppose that this is related to the strongly non-linear response of K(h) to a slight departure from full saturation. Eching and Hopmans (1993) therefore considered K s a mere fitting parameter without physical meaning in experiments that do not include full saturation. In a large soil column it is particularly difficult to establish saturation as an initial condition over the entire system. Concerning optimization quality we finally remark that our study was based on a randomized replicated design. Replication is rarely found in column experiments, probably due to cost constraints. In our case this setup allowed a proper statistical evaluation. Non-unique solutions would have increased the random error of parameter estimates due to convergence at different local minima. Statistical significance of treatment effects therefore provided an additional evidence for reliable optimization results.
Inverse estimation of hydraulic properties is based on the assumption that all relevant system properties are captured properly by the underlying hydraulic model. In our case this means that root induced changes were fully captured by changes in the optimized macroscopic Kosugi parameters (θ s , r m , σ, K s ). Implicitly we claim that other potential root effects on flow properties such as wettability (e.g. Carminati 2012; Carminati and Vetterlein 2013) and pore connectivity (e.g. Whalley et al. 2005;Feeney et al. 2006) are of minor importance for the effective hydraulic behaviour of the system.
Concerning wettability, Carminati and Vetterlein (2013) indicate a mucilage affected zone around the root of 50 μm. Underlying a maximum root length density of 4.1 cm cm −3 measured in our experiment (mustard, 0-10 cm) as well as the corresponding average root diameter of 0.265 mm, about 1.0 % of total pore volume was directly affected by mucilage effects. Here it is assumed that the entire root length had a fully functional rhizosphere sheath. This small percentage of directly affected pore space therefore suggests that changes of contact angle by mucilage were insufficient to explain the differentiation in effective flow behaviour we observed and which resulted in significantly different hydraulic parameters. Consequently we are confident that mainly structural changes of material properties, i.e. frequency and distribution of pores as expressed in the macroscopic parameters, were the predominant reason for the differences we observed between planted and unplanted soil columns. Still we recognize that structural changes and mucilage are strongly related (e.g. Tisdall and Oades 1982). Furthermore the rhizosphere affected pore space might have been higher when considering additional root length (fine roots, hair roots) which was not captured by root washing and image analysis (Pierret et al. 2005). To further elucidate the role of mucilage on different scales, both theoretically and experimentally, root architecture models and drainage experiments with fluids of different contact angle with soil (e.g. water and ethanol; Jarvis et al. 2008) might be used.
The second implicit assumption concerns pore connectivity. In our study the related macroscopic parameter (tortuosity l) of the hydraulic conductivity function was fixed. In spite of this we claimed that the key pore parameters of interest in our study (θ s , r m , σ) could be correctly estimated. It is well known that plant roots enhance pore connectivity (e.g. Pagliai and De Nobili 1993;Whalley et al. 2005). Therefore distinct estimates of tortuosity between rooted and non-rooted treatments could be considered fundamental to capture root induced changes. The main reason for not including tortuosity in the optimization was to avoid an additional free parameters and thereby the risk of non-uniqueness of the optimized solution . Tortuosity is a poorly defined fitting parameter in macroscopic models of hydraulic conductivity (Vervoort and Cattle 2003). Thus it is difficult to define proper initial values and parameter constraints. Furthermore the parameter mostly affected by an inadequate tortuosity value is K s , while our study focused on root induced changes in PSD parameters (θ s , r m , σ). Fixing the tortuosity parameter in the hydraulic conductivity function, which is potentially treatment sensitive, implies that its effect on the observed water flow is incorporated into a correlated parameter. Elliot et al. (2008) among others showed that K s is a function of pore connectivity and tortuosity. Thus mainly the optimized K s values might have been biased by an effect that would have been otherwise attributed to distinct tortuosity. The higher values of K s in the rooted columns indeed could indicate that upon optimization the potentially enhanced pore connectivity was compensated by increasing K s to properly meet the observed water flow. We recognize that parameter correlation in macroscopic hydraulic property models is a general problem for optimization (Pollacco et al. 2013). Still we consider that reduction of free parameters in the hydraulic conductivity function was more advantageous for properly estimating the water retention parameters of interest here, than obtaining a physically unclear approximation for the root effect on pore connectivity.
In this context we notice that in spite of a good statistical fit between measured and simulated outflow as shown in Table 3 it cannot be excluded that part of the physics of the system (e.g. tortuosity, contact angle) are not captured entirely. In this context the limits of macroscopic models have to be recognized. In order to account for root induced changes in pore geometry (e.g. Feeney et al. 2006) pore network models (e.g. Vogel 2000, Holtham et al. 2007) would be more appropriate. Also multi-modal models with distinct tortuosity of each domain (e.g. Dexter and Richard 2009) might be considered. Still parameterization of such models is challenging (Šimůnek et al. 2003).

Root influence on pore evolution
Modelling root driven pore dynamics first requires understanding the main changes in the PSD of rooted soil. Analysis of variance of the inversely determined hydraulic parameters revealed that there were two processes involved in the temporal change of hydraulic properties between the initial and final experiment, i.e. soil settlement and root influences. Soil settlement was expressed by the reduction of total porosity (pore closure), drift towards smaller pore classes and heterogeneization of the pore system (larger standard deviation of the PSD). These changes were more pronounced in the lower layer of the soil columns than in the upper layer. Rühle et al. (2013) recently observed a similar temporal trend towards a more heterogeneous pore system and non-uniform hydraulic behaviour during long-term soil column experiments.
Importantly we found significant differences in pore evolution between the non-planted and planted soil columns. This provided clear evidence of short term influences of roots on the PSD which have to be considered to describe the hydraulic behaviour of soil. Roots stabilized porosity, reducing pore loss (higher θ s ) due to soil settlement and enhanced heterogeneity of the pore system (increase in the pore radius standard deviations). Horn et al. (1994) and Dexter and Richard (2009) remarked that heterogenization of the pore system is an important indicator of structural porosity and results from the formation of both coarse intra-aggregate and fine inter-aggregate pores upon soil structuring. On the contrary non-planted columns showed structure degradation expressed by a narrowing of the PSD. Bodner et al. (2013) reported similar dynamics from a field study with the same treatments (bare soil, mustard, rye) where planted plots also showed a significantly higher σ.
Also an increase in K s in the rooted columns compared to the non-planted treatment was found. As shown by Hayashi et al. (2006), K s would be related to both r m as well as σ. The main reason for this is the strong influence of highly conductive and less tortuous macropores. The higher σ found in the planted columns implied that there was also an increase in the volume of high transmission macropores which strongly influenced K s (cf. Table 6). However the main volumetric change related to the root induced increase in σ was found in the micropore volume.
Having identified significant root induced changes in macroscopic parameters (θ s , σ, K s ), we could then try to identify relevant processes leading to the distinct changes in different pore classes. According to Watt et al. (2006) the volume of pores with a radius between 15 and 1,000 μm is the main space of root growth. Zobel (2008) indicated that fine roots, building up to 95 % of total root length, have been measured to a diameter as small as 60 μm. Thus we can assume that roots penetrate pores with a radius down to 30 μm, i.e. mainly the macropore range. Smaller pores can be accessed by root hairs (diameter around 10 μm; Jungk 2001) and mycorrhizal hyphae (diameter between 2 and 5 μm; Drew et al. 2003). Between 87.8 % (mustard) and 90.9 % (rye) of total root volume was within a radius range of <500 μm. Reduction in macro-and partially also mesoporosity (15 to 500 μm; 53 % of total decrease in pore volume in the rooted columns) of the rooted columns might be related to pore clogging by root ingrowth into existing pore channels. The pore space directly occupied by roots however was relatively small during the 3.5 months growing time (0.13 to 0.02 % for mustard and rye respectively). Furthermore compared to the non planted columns higher reduction in pore volume of the rooted columns was found in mesopores and larger micropores between 2.5 and 15 μm. These pores are not directly accessible to root in-growth. Therefore pore clogging can only partially explain the observed pore dynamics in actively rooted soil. Pores <2.5 μm showed a substantial increase compared to the nonplanted columns. Using the size classes given by Watt et al. (2006), roots thereby predominantly reduced the pore space between microaggregates (intra-aggregate space) while the pore space within these aggregates (inter-microaggregate porosity) increased. We suggest that capillary driven coalescence upon root water uptake (Leij et al. 2002) played a main role in the shift to higher microporosity. Also local compaction upon root penetration (Dexter 1987;Whalley et al. 2004;Aravena et al. 2011) might have increased micropore volume. Due to the dominance of fine axes in our model species however the extent of this process can be expected to be limited. The increase in pores <2.5 μm is of high functional importance due to their role in water storage. Recent results of Moradi et al. (2011) showed higher water content in the vicinity of plant roots. The authors suggested root mucilage as the main reason for this phenomenon, but also mention possible mechanical changes of soil pore configuration around roots. Although our method did not allow visualization of local phenomena, the bulk increase of storage pores in rooted soil is in agreement with this finding.
The higher value of s also resulted in an increase of larger macropores >500 μm in the rooted columns. Macropores are low in frequency but still they essentially contribute to overall porosity and hydraulic functioning of the soil. These pore class can be readily changed by localized effects in the vicinity of plant roots. The increase in macropores >500 μm in rooted columns points to such spatially heterogeneous effects. Carminati et al. (2009) described root shrinkage and air gap formation as a process of root induced macropore formation. Also fissures and cracks in the rhizosphere due to enhanced wetting and drying can act as newly formed macropores (Dexter 1987). The functional importance of these large pore classes is related the improved infiltrability and better soil aeration. The higher K s in the rooted columns however can be a response to both, volumetric increase in macroporosity as well as higher pore connectivity.
Scaling between microscopic processes at the rootpore interface and changes in effective hydraulic behaviour is still a challenge. Effective parameters of the models used here were shown to be responsive to root effects. This provides quantitative evidence of the importance of roots for hydraulic processes. The changes in different pore classes point to distinct processes of root-soil structure interaction. However macroscopic models cannot extract the relative importance of single processes involved. This might be addressed more appropriately by spatially explicit root architecture models (e.g. Leitner et al. 2010) combined to pore network models (Holtham et al. 2007).
Model framework for root induced pore evolution When studying plant roots as drivers for the development of structural porosity, a dynamic description of PSD is required. Still in hydrological modelling, soil hydraulic properties are generally considered as fixed material property that does not change over time or under the influence of any structure forming agent (Šimůnek et al. 2003). There are only few quantitative approaches that have been suggested to describe pore evolution. Or et al. (2000) and Leij et al. (2002) have published a series of papers based on a convection-dispersion like equation to model temporal changes in hydraulic properties.
This model was shown to adequately describe all cases of pore drift towards a smaller r m and a larger σ. However application of the model to the dynamics observed in our study revealed two problems. First, the model did not capture properly the dynamics in the nonplanted columns where σ showed the reverse trend. Narrowing of the pore size distribution has been described for processes of structure degradation, such as aggregate disruption or soil compaction (e.g. Starsev and McNabb 2001;Hayashi et al. 2006). While Leij et al. (2002) successfully modelled post-tillage soil settlement; tillage itself would have induced a reverse process. Thus several naturally occurring pore dynamics would require a different modelling approach. Second, our results demonstrated that the linear relation between drift and dispersion via a constant dispersivity could not capture appropriately the changes in PSD with increasing rooting density. This might be related to an inadequate drift function. The sharp drift between a nonrooted and rooted soil suggests that other mechanisms than rooting density (e.g. root exudation, non-detected fine roots or root hairs, compaction and particle rearrangement) have been involved. However our data showed that roots had a significant impact on σ while they did not significantly influence r m . In this case it was probably not adequate to functionally relate dispersion (change in σ) and drift. This indicates that there were qualitatively different processes involved or at least quantitative difference not captured by a linear relation. The pore evolution model properly characterizes a diffusion like process (shift from lower to higher entropy) equilibrating pore volume over radius with time. However this physical behaviour is questionable for the actively self-organizing biological root-microbesoil system (Young and Crawford 2004) where energy driven processes lead to a higher order in soil structure. The pore evolution model still provides a useful framework to study the nature of several root related processes. Ghezzehei and Or (2000) applied a physically based aggregate coalescence model to predict the drift and dispersion terms. Our results suggested an important role of root driven coalescence among other processes involved simultaneously over different pore ranges (clogging, compaction, stabilization). Radius dependent drift and dispersion might improve the capacity of the model to capture root effects. Still the complexity of the biologically active root-pore system makes the formulation of a comprehensive physically based model a challenge.

Conclusions
We quantified the changes in the PSD of soil rooted by mustard and rye compared to a non rooted soil. Based on a column drainage experiment parameters of the lognormal Kosugi PSD model were determined by inverse optimization. Statistical indicators showed that the inverse procedure provided reliable estimates of the macroscopic pore parameters. Although roots can influence pore connectivity, tortuosity was fixed avoid optimizing an additional fitting parameter. From the inverse results, we concluded that retention parameters can be properly estimated by macroscopic models, while the strong influence of pore geometry on the hydraulic conductivity function requires other approaches such as pore network models. We demonstrated that besides general soil settlement over time, root effects were significantly influencing temporal dynamics of the PSD. Reduction in total porosity was significantly reduced by roots. Also heterogeneity of the pore space was increased in the rooted columns indicating an increase in structural porosity. The volume of large transmission macropores as well as fine storage pore was higher in the rooted compared to the non-planted columns. From the reduction in pore space accessible to roots we concluded that pore clogging was only of minor importance, while enhanced structuring by enmeshment and aggregate coalescence were suggested as dominant processes. Applying a pore evolution model we were able to simulate the drift to smaller pore radius classes and increased pore heterogeneity in the rooted columns, while structure degradation in the non rooted columns was not captured appropriately. When using a rooting density dependent drift function however, suggested that a diffusion like model of pore evolution is only partially appropriate to describe the complex processes of active structure formation by plant roots.