Relations between macropore network characteristics and the degree of preferential solute transport

The characteristics of the soil macropore network determine the potential for fast transport of agrochemicals and contaminants through the soil. The objective of this study was to examine the relationships between macropore network characteristics, hydraulic properties and state variables and measures of preferential transport. Experiments were carried out under near-saturated conditions on undisturbed columns sampled from four agricultural topsoils of contrasting texture and structure. Macropore network characteristics were computed from 3-D X-ray tomography images of the soil pore system. Non-reactive solute transport experiments were carried out at five steady-state water flow rates from 2 to 12 mm h. The degree of preferential transport was evaluated by the normalised 5 % solute arrival time and the apparent dispersivity calculated from the resulting breakthrough curves. Near-saturated hydraulic conductivities were measured on the same samples using a tension disc infiltrometer placed on top of the columns. Results showed that many of the macropore network characteristics were intercorrelated. For example, large macroporosities were associated with larger specific macropore surface areas and better local connectivity of the macropore network. Generally, an increased flow rate resulted in earlier solute breakthrough and a shifting of the arrival of peak concentration towards smaller drained volumes. Columns with smaller macroporosities, poorer local connectivity of the macropore network and smaller near-saturated hydraulic conductivities exhibited a greater degree of preferential transport. This can be explained by the fact that, with only two exceptions, global (i.e. sample scale) continuity of the macropore network was still preserved at low macroporosities. Thus, for any given flow rate, pores of larger diameter were actively conducting solute in soils of smaller near-saturated hydraulic conductivity. This was associated with larger local transport velocities and, hence, less time for equilibration between the macropores and the surrounding matrix which made the transport more preferential. Conversely, the large specific macropore surface area and well-connected macropore networks associated with columns with large macroporosities limit the degree of preferential transport because they increase the diffusive flux between macropores and the soil matrix and they increase the near-saturated hydraulic conductivity. The normalised 5 % arrival times were most strongly correlated with the estimated hydraulic state variables (e.g. with the degree of saturation in the macropores R = 0.589), since these combine into one measure the effects of irrigation rate and the nearsaturated hydraulic conductivity function, which in turn implicitly depends on the volume, size distribution, global continuity, local connectivity and tortuosity of the macropore network.


Introduction
Preferential transport in soil occurs when water containing solutes moves predominantly through a limited part of the pore space.Preferential transport thus reduces the residence time for solutes in the unsaturated zone.The potential for preferential flow and transport is largely determined by soil properties and is generally large for structured loamy and clayey soils.In these soils, large inter-aggregate pores and biopores often act as pathways for fast transport of agrochemicals and contaminants (Thomas and Phillips, 1979).
The extent to which the potential for preferential flow and Published by Copernicus Publications on behalf of the European Geosciences Union.
transport is realised is determined by the initial and boundary conditions, where intense rainfall on initially wet soil generally leads to a high degree of preferential transport through any macropores open to the soil surface (Jarvis, 2007).Fast flow in macropores can even be generated when the soil is initially dry, especially if it has become water repellent (e.g.Jarvis et al., 2008).
The effects of flow rate on the transport of non-reactive solutes under steady-state water flow have been studied extensively using both controlled surface supply tensions (Seyfried and Rao, 1987;Langner et al., 1999;Yu et al., 2014) and controlled irrigation rates (Dyson and White, 1989;Haws et al., 2004;Pot et al., 2005).Generally the resulting solute breakthrough curves are shifted towards a smaller number of drained pore volumes and more skewed with higher flow rates and greater saturation indicating a higher degree of preferential transport.However, for some soils the water flow rate seems to have little or no effect on the shape of solute breakthrough curves (Schulin et al., 1987;Dyson and White, 1989).All of the above cited studies have used models to evaluate the degree of preferential transport.A common approach has been to fit both the convectiondispersion equation (CDE) and the mobile-immobile model (MIM) to the breakthrough curves (Seyfried and Rao, 1987;Langner et al., 1999;Yu et al., 2014).Equilibrium transport has been assumed if the fit of the CDE to the measured breakthrough curve is as good as the MIM.Optimised values for the parameters describing the immobile water content and the mass transfer rate between the mobile and immobile domains have been used as indicators of preferential transport (Langner et al., 1999).This approach is justified by the fact that the CDE is derived from assumptions of perfect lateral mixing (i.e.zero preferential transport).However, it has been shown that the CDE is very flexible and, hence, can be well fitted to breakthrough curves with a very fast breakthrough (Koestel et al., 2011).Deviations from the CDE may, therefore, not be a suitable indicator of preferential transport from a leaching risk perspective.Furthermore, the parameters in both the CDE and the MIM are often highly correlated which means that many parameter combinations may give equally good fits to measured data (Koch and Flühler, 1993;Beven et al., 1993).Physical interpretations of fitted parameter values are, therefore, difficult.An attractive alternative to this approach is to evaluate the degree of preferential transport using model independent shape measures or indicators (Kamra and Lennartz, 2005;Koestel et al., 2011).
In many studies to date, indirect indicators of soil structure observed or measured at larger scales have been related to solute breakthrough curves (Vervoort et al., 1999;Jarvis, 2007) due to the lack of efficient experimental technologies to directly quantify structure at the pore scale.One method to quantify macropore networks is through analysis of 2-D sections of resin-impregnated soil.However, this method is labour intensive and the number and size of investigated samples have, therefore, been limited.Holland et al. (2007) used this method to calculate pore network characteristics from 2-D images of horizontal soil sections and related these measures to breakthrough curves from nonreactive solute transport experiments carried out on separate columns.They fitted log-normal transit time distributions to the breakthrough curves and showed that breakthrough curves for columns with a less well connected pore network had smaller mean transit times and larger standard deviations of transit times which we interpret as a higher degree of preferential transport.These columns also had smaller macroporosities and specific macropore surface areas.In recent years, non-destructive analysis of soil macropore networks has become popular through different 3-D imaging techniques with X-ray tomography being the most common (Wildenschild and Sheppard, 2013).Quantitative measures of the macropore network (e.g.macroporosity, connectivity and tortuosity) calculated from X-ray tomography images can provide useful information on soil structure, although they are dependent on the original image quality, image processing methods and the algorithms used (Schlüter et al., 2014).Thus, X-ray tomography has the potential to link soil structural features to soil functioning.Although only a few studies have been reported, X-ray data has led to some new insights in the study of hydraulic and transport processes in soil.For example Vanderborght et al. (2002) used X-ray tomography to characterise the macropore network of a loamy forest soil.They showed that a dense macropore network at the 10-20 cm depth resulted in homogeneous transport whereas isolated large continuous macropores which were present at a depth of 50-60 cm resulted in more heterogeneous transport.Luo et al. (2010) showed that macropore network characteristics such as the hydraulic radius, the number of continuous flow paths and flow path angle were correlated to dispersivities calculated from fitting the CDE to tracer breakthrough curves obtained under saturated conditions.They explained large dispersivity values, indicative of a high degree of preferential transport, with the presence of only one or a few large continuous macropores.Soil columns with larger and better connected macropore networks showed slower solute breakthrough.
For many soils and climate conditions much of the transport takes place when the topsoil is not fully saturated.As far as we know, the relationships between macropore network characteristics derived from 3-D images of the pore system and solute transport under unsaturated conditions at a range of different flow rates have not yet been studied.The objective of this study was to determine how macropore network characteristics influence the degree of preferential transport of a non-reactive tracer under near-saturated conditions in agricultural topsoils.This was achieved by combining X-ray tomography with measurements of near-saturated hydraulic conductivity using tension disc infiltrometers and non-reactive solute transport experiments carried out at different steady-state water flow rates in four soils of contrasting texture and structure.

Soils
Soil columns were sampled in PVC pipes (20 cm high, 20 cm diameter) using a tractor mounted hydraulic press at four conventionally tilled fields close to Uppsala in eastern Sweden.Five replicate samples were taken from the topsoil at each site leaving the soil surface undisturbed.The soils were Säby 1 (loam), Säby 2 (clay), Krusenberg (clay loam) and Ultuna (clay) (Table 1).These soils were chosen because they have clay contents above 20 % and, hence, were expected to have a well-developed structure (Horn et al., 1994) and a potential for preferential transport (Koestel and Jorda, 2014).The pore systems of cultivated topsoils change with time, for example due to tillage and subsequent consolidation and swelling/shrinking.In this study, all samples were taken on 26 June 2013.Three of the sampling sites (Säby 2, Krusenberg and Ultuna) were sown with spring crops a few weeks earlier (Table 1).At the time of sampling, the seedbeds for these soils consisted of fine individual aggregates.In contrast, desiccation cracks extending a few cm down into the soil were visible for the Säby 1 soil which had been sown with an autumn crop.The samples were stored at 3 • C until the start of the solute breakthrough experiments.Since the natural soil surface was uneven, the length of the soil columns varied between 15.0 and 18.3 cm assuming that the soil surface was located where the soil matrix, as determined by X-ray tomography, filled 50 % of the horizontal column cross-section area.
One of the Säby 2 (clay) columns contained a large amount of straw residues and was accidentally destroyed during transport and one of the Ultuna (clay) columns ponded already at the lowest irrigation rate.These columns were excluded from the study.Hence, transport experiments, X-ray tomography and measurements of near-saturated hydraulic conductivity were carried out on 18 out of the 20 columns (5 columns for Säby 1 (loam) and Krusenberg (clay loam) and 4 columns for Säby 2 (clay) and Ultuna (clay)).

Non-reactive solute transport experiments
Solute breakthrough experiments were conducted in an irrigation chamber at the Department of Soil and Environment at the Swedish University of Agricultural Sciences (SLU), during October and December 2013.The irrigation chamber is described in detail in Liu et al. (2012).Before the start of the experiments polyamide cloths (mesh size 50 µm) were attached to the bottom of the samples in order to minimise soil loss.The irrigation nozzles were turned on each minute and the flow rates were adjusted by setting the time each individual irrigation nozzle was on during 1 minute.To achieve the desired irrigation rates calibration was carried out before the start of each new irrigation.During the experiments the effluent water volumes were monitored daily.Irrigation in-tensities were not always constant at the nominal rates possibly due to the formation or dissolution of precipitates in the nozzles during the experiment.The actual average irrigation rates and their standard deviations estimated from effluent volumes are presented in Table S1 in the Supplement.
The columns were first irrigated for 6 days with tap water at an intensity of 2 mm h −1 prior to the first solute application in order to allow the soil to settle and to achieve stable flow rates and background effluent electrical conductivities.Five tracer experiments were carried out under subsequently increasing steady-state irrigation rates of 2, 4, 6, 8 and 12 mm h −1 .As irrigation intensities increased the saturated hydraulic conductivities of a few columns were exceeded and surface ponding occurred.When a column ponded the transport was not further analysed.After an increase in irrigation rate, a warm-up period was needed for the flow rates and background conductivities in the column effluents to stabilise.When the system was at approximate steady-state a 2 ml pulse of potassium bromide solution (250 mg bromide added per ml of tap water) was applied during approximately 30 s at the surface of each column using a pipette.A reasonably homogeneous spatial application was achieved using a plastic grid placed on the surface before application.For all irrigation rates at least 185 mm water had passed through the columns when the experiments were stopped.
The electrical conductivity of the effluent was measured in flow-through vessels (D201, WTW GmbH, Weilheim, Germany) using electrical conductivity meters (Cond 3310, WTW GmbH, Weilheim, Germany).Measurements were recorded at 5 min resolution for the irrigation intensities 2-8 mm h −1 and at 1 min resolution for the 12 mm h −1 irrigation rate.Background electrical conductivities in the water supply were measured in grab samples taken approximately every 72 h throughout the experiment from one of the irrigation nozzles.The background electrical conductivities varied between 473 µS cm −1 and 516 µS cm −1 .Furthermore, the background electrical conductivity of the effluent solution was influenced by resident ions in soil solution that reacted with the water flowing through the soil.We estimated the background concentration in the effluent solution during the experiments, EC background , by fitting an equation of the following form: where p 1 to p 4 are fitted parameters, to the measured electrical conductivities during the time period prior to bromide application.This equation was chosen because it gave an acceptable fit to measured data (average coefficient of determination = 0.78).Bromide concentrations in the effluent water, C out (mg ml −1 ), were calculated from the electrical conductivity data assuming a linear relationship between bromide concentration and electrical conductivity (CRC Handbook of Chem- istry and Physics, 1989): where C in (mg ml −1 ) is the known bromide concentration in the applied bromide solution, EC in (µS cm −1 ) and EC out (µS cm −1 ) are the electrical conductivities in the applied solution and in the effluent water respectively and t i are the times when the electrical conductivity was recorded.
After the solute transport experiments the soil columns were left to drain for 1 week and then stored at 3 • C until the X-ray tomography was carried out.We assumed that changes in the macropore system were negligible during this period of drainage.

Measures of preferential transport
Preferential transport is generally associated with an early arrival of solutes and pronounced tailing of solute breakthrough curves (Brusseau and Rao, 1990).In this study, we used two shape measures to indicate the degree of preferential transport: the normalised 5 % arrival time which reflects the tendency for early arrival of the solutes and the Eulerian apparent dispersivity, hereafter referred to as the apparent dispersivity, which is a measure of the spread of arrival times (Knudby and Carrera, 2005;Koestel et al., 2011).The normalised 5 % arrival time has previously been shown to be strongly correlated to other model independent measures that reflect early solute arrival (e.g. the holdback factor and the arrival time for the maximum concentration; Koestel et al., 2011;Ghafoor et al., 2013).
Traditionally, solute breakthrough curves have been analysed in relation to the number of pore volumes that have passed through the sample.However, for unsaturated conditions the volumetric water content, θ (−), and, hence, the pore volume are usually unknown.We, therefore, used the solute breakthrough curve to define a "tracer-specific water filled pore volume" θ eff as follows: where q (mm h −1 ) is the water flux per column cross-section area, L (mm) is the length of the soil column and µ 1 (h) is the normalised first temporal moment of the breakthrough curve (i.e. the average arrival time of the solute) given by where t (h) is time and m 0 (−) and m 1 (h) are the zeroth and first temporal moments, respectively.The integrals in Eq. ( 4) were estimated directly from the raw data.A normalised breakthrough curve can be derived as follows (Jury and Roth, 1990): where T = t µ 1 (−) is the number of effective pore volumes that has passed through the sample.The normalised cumulative breakthrough curve is defined by the following: The normalised 5 % arrival time, p 0.05 , was finally calculated from the following: The 5 % arrival time is illustrated in Fig. 1.The apparent dispersivity which is a frequently used measure of the spread in solute arrival times (e.g.Vanderborght et al., 2001;Koestel et al., 2011) was calculated according to the following: where µ 2 is the second central moment of the breakthrough curve given by the following: To enable a proper comparison between the breakthrough curves they were truncated at a time corresponding to 185 mm of applied water, which was the smallest amount that was applied in the experiments.Truncated breakthrough curves will generally result in overestimations of p 0.05 and underestimations of λ app compared to modelled breakthrough curves which allow for extrapolation beyond the measured data (Koestel et al., 2011).

X-ray computed tomography
In this study, we used the GE Phoenix v|tome|x m X-ray tomograph installed at the Department of Soil and Environment at SLU.It has a 240 kV X-ray tube, a tungsten target (beryllium window) and a GE 16 flat panel detector.We collected 2000 radiographs per column with a discretization of 2024 × 2024 pixels, corresponding to a resolution of 121.1 µm.The X-ray scans were carried out at a voltage of 200 kV with a current in the range 350-420 µA varying between columns.The exposure time for each radiograph was 250 µs.The radiographs were subsequently inverted to a 3-D image using the GE image reconstruction software datos|x and exported as TIFF-stacks (tagged image file format) with 16-bit greyscale resolution.The resulting spatial resolution of the reconstructed 3-D images was 121.1 µm in all directions.The X-ray tomography was carried out during December 2013 and February 2014, after the solute transport experiments.

Image processing and analysis
Image processing and analysis were carried out using the Fiji distribution (Schindelin et al., 2012) of the software ImageJ (Abramoff et al., 2004) including the plugin BoneJ (Doube et al., 2010), GeoDict (http://www.geodict.com)and R (R Core Team, 2012).First, the image resolution of each frame was reduced by a factor of two in all dimensions to limit the computation time during the following image processing steps.As a result, each image voxel had an edge length of 242.2 µm.

Definition of the sample volume
In order to calculate the soil macroporosity the total volume of the soil inside the PVC cylinder (hereafter referred to as the sample) needed to be estimated.Since the samples were taken from the soil surface their volumes were not identical.In order to define the sample volume, images were first rotated in order to get the columns in a perfectly upright position.Then a circular region of interest corresponding to the inner diameter of the soil columns was defined and all image sections not containing soil were cut away.The 3-D image was imported into R as a sequence of vertical slices of the original image.The soil surface was determined by searching the first increase in grey value in each XY coordinate that exceeded 25 % of the total range of grey values at that XY coordinate, starting from the top of the soil column.All voxels above and below the soil surface were then set to 0 and 1, respectively.Finally we applied a 2-D mean filter with a radius of 6 pixels to all zero valued pixels in each horizontal slice.Values below 0.7 were then set to 0 while values above 0.7 were set to 1.This filtering filled surface vented pores with a radius smaller than approximately 2 mm whilst preserving the edges of larger surface vented pores.The result of the above described operations was a binary image containing the sample volume.

Illumination correction and segmentation
The large column diameters and large soil bulk densities resulted in differences in illumination between the centre of the column and the perimeter.In addition, ring artefacts were clearly visible at approximately 15 cm depth.These artefacts, caused by beam hardening, were reduced by first assigning all voxels in the sample volume into categories based on their vertical coordinate, z, and radial coordinate, r, (i.e.their distance from the centre of the column in the horizontal direction).All grey values, GV x,y,z , were then corrected for radial differences in illumination according to the following: where GV corr,x,y,z are the corrected grey values, Perc 60,ref and Perc 80,ref are reference values for the 60th and 80th percentile grey values, respectively and Perc 60,r,z and Perc 80,r,z are the 60th and 80th percentiles of GV x,y,z at depth z and radial distance r.These percentiles were chosen because it was clear from a visual inspection of the 3-D images that the sum of the macroporosity and the fraction of organic matter at any given radius from the centre of the column in a horizontal slice never exceeded 60 % and that the stone content never exceeded 20 %.It could, therefore, be assumed that Perc 60,r,z and Perc 80,r,z represented the soil matrix and were only marginally influenced by differences in macroporosity (Iassanov and Tuller, 2010).The reference values for the 60th and 80th percentiles were chosen arbitrarily as the measured 60th and 80th percentiles for one of the columns.
The reference values only determine the absolute scale of the grey values and have no consequences for the subsequent binarization which is independent of the absolute scale.Illumination corrections were done in R. The corrected images were filtered using the 3-D median filter with a radius of 2 in ImageJ to reduce noise.A consequence of this filtering is the removal of the smallest pores.Finally, the images were binarised into macropores and soil matrix using Otsu's thresholding method (Otsu, 1979) applied to each slice.Hence, we defined macropores operationally as all pores visible in the binary images.Since the image resolution was 242.2 µm, pores with a diameter larger than approximately 484.4 µm were visible.

Measures of the macropore network
In addition to the total pore volume, a complete geometric and topological description of the pore network in soil should include metrics accounting for the size distribution of pores, their specific surface area and local connectivity (Vogel et al., 2010) as well as their heterogeneity and global (i.e.samplescale) continuity (Renard and Allard, 2013).From the binary images of the pore networks we calculated a number of such measures of both the total macropore network and the largest macropore cluster connecting the soil surface to the bottom of the sample (Table 2).The macropore clusters were determined using the "Analyse particles" ImageJ plugin, which is an adaptation of the 3-D object counter described in Bolte and Cordelières (2006).To estimate macropore sizes we calculated the pore "thickness" using BoneJ.The pore thickness is defined for each macropore voxel as the diameter of the largest sphere that fits into the macropore and contains the voxel.The histograms of thickness values were used as estimates of pore size distributions.Thicknesses were also calculated for the soil solid matrix as a proxy for aggregate width.
The existence of a pore cluster connecting the soil surface to the bottom of the sample was used as an indicator of global continuity at the scale of the sample, while the Euler number for the largest connected cluster was used as an estimate of local pore connectivity (Vogel and Kretzschmar, 1996).The mass fractal dimension was calculated as a measure of the heterogeneity of the spatial distribution of macroporosity (Peyton et al., 1994;Perret et al., 2003).The smallest pore neck on the critical path connecting the soil surface to the bottom (the critical pore diameter) and the length of this path were estimated using the PoroDict module in GeoDict.The critical macropore diameter is estimated as the diameter of the largest sphere that can pass through the pore space from the top to the bottom of the column.

Unsaturated hydraulic conductivity
After X-ray scanning, we measured the steady-state infiltration rate at the soil surface of each column with a tension disc infiltrometer set at supply tensions of ψ =5 cm and ψ =1 cm.The bottom plate of the tension infiltrometer had a diameter of 15 cm leaving space for air to escape through the soil surface near the column walls.A layer of moist fine sand was applied at the soil surface to ensure good contact between the disc and the soil.We used the steady-state infiltration rate as an estimate of the unsaturated hydraulic conductivity, K (mm h −1 ), at the supply tension given the 1dimensional nature of the flow system.The steady-state pressure potentials attained in the soil at the fluxes applied during the solute transport experiments were calculated assuming a linear relationship between log(ψ) and log(K) (Fig. 2).Such a relationship has been fitted with a coefficient of determination > 0.9 for 90 % of the data entries included in a global database containing 753 individual data sets (Jarvis et al., 2013).From the steadystate pressure potential we calculated the size of the largest water-filled pore during the solute transport experiments using the Young-Laplace capillarity equation.We then estimated the degree of saturation in the macropores at each flow rate from the X-ray derived pore size distribution.

Statistics
Spearman rank correlation coefficients were calculated between macropore characteristics, hydraulic conductivities and measures of preferential transport.Selected correlations were also graphically displayed and analysed by linear regression.For both cases, results were considered significant if p values were less than 0.05.

Soil characteristics
Figure 3 shows example images of the macropore networks generated by X-ray tomography for the four soils.All soils had larger macroporosities in the top 5 cm which corresponds to the harrowed layer.Macroporosities estimated from X-ray tomography images were between 3.3 and 12 % (Fig. 4) with an average value of 6.5 %.These values are large compared to macroporosities derived from X-ray tomography images reported in the literature for agricultural topsoils (Kim et al., 2010;Luo et al., 2010;Naveed et al., 2013).It should be noted that three of the soils included in this study were sampled a few weeks after ploughing and seedbed preparation, operations which generally increase macroporosity.Near saturated hydraulic conductivities at tensions of 1 and 5 cm, K 1 and K 5 , were in the ranges 0.35 to 38 mm h −1 and 0.14 to 0.72 mm h −1 , respectively (Fig. 4), which is in line with previously reported values for loamy and clayey agricultural topsoils (Jarvis et al., 2013).Macropore network characteristics are presented in Tables S2 and S3 (Supplement).Figure 5 shows that many of the macropore network characteristics were inter-correlated.Here we briefly discuss the implications of some of these correlations.Figure 5 shows that total macroporosity was significantly correlated with K 1 but not with K 5 , which is mainly determined by pores that are not visible at the image resolution.Luo et al. (2010) and Kim et al. (2010) also found significant correlations between saturated hydraulic conductivity and X-ray imaged macroporosity.The positive correlation between macroporosity and fractal dimension shows that the macropore networks in columns with large macroporosities are spatially more homogeneously arranged (i.e. they "fill" the total volume to a larger extent).This spacefilling property of columns with large macroporosities is also indicated by the strong negative correlation between macroporosity and mean aggregate width.Dal Ferro et al. (2013) reported positive correlations between macroporosity and fractal dimension for a silty to sandy loam soil under different fertilisation regimes.Positive correlations have also been reported between depth distributions of macroporosity and the  fractal dimension for a sandy loam soil (Perret et al., 2003) and a silty loam soil (Luo and Lin, 2009).
A positive correlation was found between macroporosity and the fraction of the total pore space consisting of the largest pore cluster.As might be expected, the Euler number, as a measure of local connectivity, showed a significant negative correlation with the macroporosity of the largest pore cluster and the fraction of the total macroporosity which was found in this cluster, although the negative correlation between the Euler number and total macroporosity was not significant.These results indicate that large macroporosities are associated with a more connected macropore network.This relationship is also suggested by the global connectivity measure, C, which indicates whether any pore cluster extends from the top to the bottom of the column (C = 1) or not (C = 0).The two columns where such a cluster did not exist had the second and third smallest macroporosities out of the 18 columns (Table S2, Supplement).

Solute transport
Figure 6 shows that, as reported in previous studies (Koestel et al., 2012;Ghafoor et al., 2013), the two measures describing the shape of the breakthrough curve, the normalised 5 % arrival time and the apparent dispersivity, were strongly correlated to each other, especially for columns exhibiting only weak preferential transport.Therefore, in the following we mainly discuss the normalised 5 % arrival times, which generally showed slightly stronger correlations with the macropore network characteristics, state variables and hydraulic properties of the soils.
All solute breakthrough curves are presented in Figs.S1 to S4 (Supplement).Generally, an increased flow rate resulted in earlier breakthrough and a shifting of the peak concentration towards smaller drained volumes which is in line with previous studies (e.g.Haws et al., 2004;Pot et al., 2005;Koestel et al., 2012).The relationship between peak height and the flow rate was less clear (Figs.S1-S4, Supplement).We are not aware of any studies that discuss how the peak height is related to macropore network characteristics.Although peak heights may contain interesting information regarding the transport process they were not further considered in this study.
The normalised 5 % arrival times varied between 0.04 and 0.52 for the four soils (Table S4 in the Supplement).Fortythree percent of the breakthrough curves had normalised 5 % arrival times smaller than the median value of approximately 0.25 presented in Koestel et al. (2012) for 302 breakthrough curves from undisturbed soil columns sampled from arable land.Thirty-nine percent of the breakthrough curves had apparent dispersivities larger than the median value of approximately 10 cm for the same data set.The decreasing number of data points with increasing flow rate for Säby 1 (loam), Säby 2 (clay) and Ultuna (clay) is due to surface ponding.The pattern seen in Figs.S1-S4 (Supplement), with decreasing normalised 5 % arrival times at higher flow rates, was not apparent in the average values for the two clay soils.This is because the columns that ponded had smaller normalised 5 % arrival times than the average values.Figure 7 shows that it is not possible to make a general statement comparing the degree of preferential transport in different soils based on breakthrough curves from one flow rate only.For exam-Figure 5. Spearman (rank) correlation coefficient between pore system characteristics, state variables and measures of preferential transport where ε mac,tot , total macroporosity; SA, specific macropore surface area; r hyd , hydraulic radius; Th mean,tot,p , mean thickness of macropores; Th mean,tot,a , mean thickness of soil matrix; D, fractal dimension; C, global connectivity; f clust , fraction of cluster volume of total macropore volume; Th mean,clust,p , mean thickness of macropores for cluster; E, Euler number; d critical , critical pore diameter; ξ , tortuosity; K 1 and K 5 , hydraulic conductivities at 1 and 5 cm tensions, respectively; 5 % Xmm h −1 , 5 % arrival times at flow rates X mm h −1 .Stars indicate significant correlations (p value < 0.05).ple at an irrigation rate of 2 mm h −1 Säby 2 (clay) exhibits stronger preferential transport than Säby 1 (loam), although this is not statistically significant, while the opposite is true for irrigation rates higher than 4 mm h −1 .

Relationships between macropore network characteristics and measures of preferential transport
Figure 8 shows the relationships between pore network characteristics and measures of preferential transport.Significant positive correlations were found between the normalised 5 % arrival times and macroporosity at irrigation rates between 2 and 6 mm h −1 (Fig. 8a).These results indicate that preferential solute transport was significantly weaker in the columns with large macroporosities, which also contained more homogeneously distributed and better-connected macropore networks and had larger hydraulic conductivities close to saturation (Fig. 5).At first sight, these results may seem counter-intuitive.However, this relationship can be explained by the fact that macropores will fill up until the hydraulic conductivity of the soil is large enough to conduct water at the set flow rate.For soils with small macroporosities a larger fraction of the macroporosity must be water filled (i.e. the degree of saturation is larger) at the given flow rate (assuming similar pore size distributions).This means that pores with larger diameter, which conduct water at higher flow rates, will be active in the transport.With less time for equilibration of solute concentrations between the macropores and the surrounding matrix the transport becomes more preferential.The important role of diffusion in controlling preferential transport in macroporous soils has previously been discussed in a number of studies (e.g.Jarvis et al., 1991;Gerke and van Genuchten, 1993;Koestel and Larsbo, 2014).
In our study, normalised 5 % arrival times were positively correlated to the specific macropore surface area and negatively correlated to mean aggregate thickness (Figs. 5 and  8c).These two macropore network characteristics were in turn negatively correlated to each other.For an ideal aggregate geometry (e.g.uniform sized cubes or spheres), they would be perfectly correlated.These results suggest that a larger area for diffusive exchange between macropores and the soil matrix increases the diffusive flux, thereby limiting preferential transport.The significance of aggregate thickness may also partly reflect the effects of pore connectivity on preferential transport, since soils with smaller aggregates also have a denser and better connected macropore network compared to soils with larger aggregates (Fig. 5).
The Euler number was negatively correlated to the normalised 5 % arrival times, although only significantly for the two lowest irrigation rates (Fig. 8g).Thus, it seems that provided global continuity is maintained, which was the case for all but two columns in this study, a less connected pore network leads to stronger preferential transport.These results are in line with previous studies (Holland et al., 2007;Luo et al., 2010) and support the conceptual model described by Jarvis (2007) and further elaborated on by Jarvis et al. (2012), where the strongest preferential transport was predicted to occur in soils which contain large continuous macropores, but lack well-connected networks of smaller macropores and mesopores.A well-connected network of smaller pores reduces the likelihood of soil water potentials increasing sufficiently to trigger flow in the largest macropores and increases rates of lateral convective and diffusive mass exchange during vertical flow.At first glance, these findings may seem to contradict the results of numerical studies carried out at the macroscopic Darcy scale, where it has been shown that fields containing one or a few paths of high hydraulic conductivity connected in the flow direction exhibit stronger preferential flow than fields with poorly connected zones of high hydraulic conductivity (e.g.Knudby and Carrera, 2005;Bianchi et al., 2011).However, there is no conflict between these findings at the pore and Darcy scales.This is illustrated schematically in Fig. 9, which shows that provided pore continuity across the sample is maintained, decreased macroporosity and associated poor local connectivity at the pore scale should result in stronger preferential transport.At the macro- scopic Darcy scale, this will be manifested as spatially correlated high conductivity flow pathways surrounded by low conductivity zones.The conceptual model presented in Fig. 9 can also be interpreted in a percolation theory framework (Berkowitz and Balberg, 1993;Hunt, 2001).Preferential flow reaches a maximum as we approach the percolation thresh-old from above (i.e.moving from Fig. 9a to b), while Fig. 9c illustrates a discontinuous macropore network below the percolation threshold that cannot sustain preferential flow.
The correlations between pore network characteristics and measures of preferential transport were generally stronger at the lower irrigation rates (Figs. 5 and 8).One reason for this Figure 9. Schematic illustration of the relationship between local connectivity of the pore network, macroscopic (Darcy scale) hydraulic conductivity and preferential transport (a) this represents a well-connected macropore network (a single pore cluster) with a uniformly high macroscopic hydraulic conductivity denoted by dark shading.The degree of preferential transport is consequently low (b) this network was created from the a-network by removing some pores to create seven isolated pore clusters.It represents a macropore network which is poorly connected locally, but which exhibits a high degree of preferential transport because one of the pore clusters is still continuous through the sample.At the Darcy scale, this is reflected in a continuous high conductivity pathway through the sample (dark shading) surrounded by low conductivity matrix (light shading) (c) this shows that removing one more pore results in a breakdown of sample-scale continuity, which will limit preferential transport.Macroscopic hydraulic conductivity is now uniformly low.
is simply that there were fewer columns for the higher irrigation rates due to surface ponding.Another tentative explanation is that the calculated pore network characteristics depend mostly on the smaller macropores that dominate transport at lower flow rates.For example the Euler number is mainly determined by the dense network of finer pores while the transport at higher irrigation rates is to a larger extent determined by the larger macropores which comprise a smaller fraction of the total imaged pore volume.

Relationships between hydraulic properties and state variables and measures of preferential transport
Figure 10 shows the relationships between hydraulic properties and state variables and measures of preferential transport.The positive correlation between K 1 and 5 % arrival times was significant for all irrigation rates (Fig. 10a).This shows that columns with large near-saturated conductivities exhibit a low degree of preferential transport.This was also found by Ghafoor et al. (2013) for a large data set of 45 topsoil columns sampled from 14 agricultural fields in a small catchment in south-east Sweden.Again, this relationship arises because at any given applied irrigation rate, larger pores will become water-filled in soils of lower near-saturated hydraulic conductivity.This will tend to enhance the degree of preferential transport providing macropore continuity does not break down at the sample scale.This is further illustrated by Fig. 10c and d, which shows relationships between the degree of preferential transport and the estimated size of the largest water-filled macropore (p value < 0.001, coefficient of determination = 0.485 and 0.438 for all data in Fig. 10c and d, respectively).A similar relationship was previously reported by Ghafoor et al. (2013).Since the near-saturated hydraulic conductivity is dependent on the pore volume as well as the continuity, connectivity and tortuosity of the macropore network these properties are implicitly accounted for in the calculations of the largest water-filled pores.Plotting the data in this way greatly reduces the scatter due to variations in irrigation rate, which is apparent in the plots of breakthrough shape measures against pore network characteristics (Fig. 8) and hydraulic conductivity (Fig. 10a and b).This is even more apparent in the relationships between the degree of saturation in the macropores and the normalised 5 % arrival times and apparent dispersivity shown in Fig. 10e and f (coefficient of determination = 0.589 and 0.520 for all data in Fig. 10e and f, respectively).In addition to the diameter of the largest water-filled macropores, it is reasonable to assume that the fraction of macropores with a diameter close to the diameter of the largest water-filled pore in relation to the total waterfilled macropore volume (i.e. the pore size distribution of the water-filled macropores) should affect the degree of preferential transport.A strong correlation between the degree of saturation in the macropores and normalised 5 % solute arrival times was also found by Koestel et al. (2013) for 65 topsoil columns sampled from one agricultural field in Denmark.
There are many possible reasons for the remaining unexplained scatter in the data.Most importantly the transport behaviour of these soils should not be expected to be fully determined by the degree of saturation in the macropores.As discussed above the specific macropore surface area will affect the degree of preferential transport through its influence on the diffusive exchange between pore domains (Fig. 8c and  d).Characteristics of the pores with sizes below the resolution of the images are also likely to influence solute transport.The estimated size of the largest water-filled pores and, hence, the degree of saturation is dependent on the validity of the assumption of a linear relationship between log(ψ) and log(K).The calculated macropore characteristics are estimates of the true soil properties whose accuracy depends on the quality of the images and the image processing procedure.The estimated normalised 5 % arrival times and apparent dispersivities are dependent on the accuracy of the subtraction of background electrical conductivity values and were possibly also influenced by fluctuations in the irrigation rate (Table S1 in the Supplement).Finally, it can be noted that the macropore networks are not constant with depth (Fig. 3) and the use of averages for the whole sample may, therefore, not be appropriate.Considering all these uncertainties and potential sources of error we consider that the relationships shown in Figs. 8 and 10 are remarkably strong.water-filled pore, the degree of saturation in the macropores and measures of preferential transport 5 for the four soils.Linear regression lines (dashed) were drawn if the slope was significantly (p-6 value<0.05)different from zero.7 8 Figure 10.Relationship between hydraulic conductivity at ψ = 1 cm, the diameter of the largest water-filled pore, the degree of saturation in the macropores and measures of preferential transport for the four soils.Linear regression lines (dashed) were drawn if the slope was significantly (p value < 0.05) different from zero.

Conclusions
The results of this study indicate that soils with large macroporosities are associated with a well-connected, space-filling macropore network with a large specific macropore surface area.These properties limit the degree of preferential transport because they increase the diffusive flux between macropores and the soil matrix and increase the near-saturated hydraulic conductivity, thereby inactivating the largest pores during solute transport under steady near-saturated flow.Our results suggest that the near-saturated hydraulic conductivity function combined with some measure reflecting infiltration rates (i.e. the upper boundary condition) could be used to estimate preferential transport under field conditions, because together they determine the size of the largest water-filled pore.This study was carried out on conventionally tilled agri-cultural topsoils with clay contents above 20 % using constant irrigation intensities from 2 to 12 mm h −1 which resulted in nearly saturated conditions.The extent to which the conclusions of our study with respect to the relationships between pore network characteristics and flow and transport properties can be extrapolated to other soil and climatic conditions should be the subject of further research.However, it is encouraging that our results are in agreement with the general conceptual model proposed by Jarvis et al. (2012) and are also consistent with the results of the only three previous studies we are aware of that have investigated these relationships (Vanderborght et al., 2002;Holland et al., 2007;Luo et al., 2010).

Figure 1 .
Figure 1.Illustration of the normalised 5 % arrival time, p 0.05 .(a) shows the travel-time probability density function, g n , where the area to the left of the 5 % arrival time equals 5 % of the total area under the curve and (b) shows the corresponding cumulative distribution function of the travel times, G n .Both curves are plotted against the number of effective pore volumes drained (T ).

Figure 2 .
Figure 2. Determination of soil tensions corresponding to steadystate irrigation rates of 2 and 12 mm h −1 for two of the columns

M.
Larsbo et al.: Macropore network characteristics and the degree of preferential solute transport

Figure 4 .
Figure 4. Macroporosities and near-saturated hydraulic conductivities for the four soils.

Figure 6 .
Figure 6.Correlations between measures of preferential transport.

Figure 7 .
Figure 7. Mean 5 % arrival times for the four soils.The number of breakthrough curves, n, was smaller at higher irrigation rates due to ponding.Error bars indicate 1 standard deviation

Figure 8 .
Figure 8. Relationship between macropore characteristics and measures of preferential transport for the four soils.Linear regression lines (dashed) were drawn if the slope was significantly (p value < 0.05) different from zero.

Figure 10 .
Figure 10.Relationship between hydraulic conductivity at ψ = 1 cm, the diameter of the largest 4water-filled pore, the degree of saturation in the macropores and measures of preferential transport 5 for the four soils.Linear regression lines (dashed) were drawn if the slope was significantly (p-6 value<0.05)different from zero.7

Table 2 .
Measures of the macropore network.i.e. the shortest path from the surface to the bottom with pore diameters larger than or equal to the critical pore diameter) divided by the length of the sample