Prediction of biopore-and matrix-dominated flow from X-ray CT-derived macropore network characteristics

Prediction and modeling of localized flow processes in macropores is of crucial importance for sustaining both soil and water quality. However, currently there are no reliable means to predict preferential flow due to its inherently large spatial variability. The aim of this study was to investigate the predictive performance of previously developed empirical models for both water and air flow and to explore the potential applicability of X-ray computed tomography (CT)-derived macropore network characteristics. For this purpose, 65 cylindrical soil columns (6 cm diameter and 3.5 cm height) were extracted from the topsoil (5 cm to 8.5 cm depth) in a 15 m× 15 m grid from an agricultural field located in Silstrup, Denmark. All soil columns were scanned with an industrial X-ray CT scanner (129 μm resolution) and later employed for measurement of saturated hydraulic conductivity, air permeability at −30 and −100 cm matric potential, and gas diffusivity at −30 and −100 cm matric potential. Distribution maps for saturated hydraulic conductivity, air permeability, and gas diffusivity reflected no autocorrelation irrespective of soil texture and organic matter content. Existing empirical predictive models for saturated hydraulic conductivity and air permeability showed poor performance, as they were not able to realistically capture macropore flow. The tested empirical model for gas diffusivity predicted measurements at −100 cm matric potential reasonably well, but failed at −30 cm matric potential, particularly for soil columns with biopore-dominated flow. X-ray CT-derived macroporosity matched the measured airfilled porosity at −30 cm matric potential well. Many of the CT-derived macropore network characteristics were strongly interrelated. Most of the macropore network characteristics were also significantly correlated with saturated hydraulic conductivity, air permeability, and gas diffusivity. The predictive Ahuja et al. (1984) model for saturated hydraulic conductivity, air permeability, and gas diffusivity performed reasonably well when parameterized with novel, X-ray CTderived parameters such as effective percolating macroporosity for biopore-dominated flow and total macroporosity for matrix-dominated flow. The obtained results further indicate that it is crucially important to discern between matrixdominated and biopore-dominated flow for accurate prediction of macropore flow from X-ray CT-derived macropore network characteristics. Published by Copernicus Publications on behalf of the European Geosciences Union. 4018 M. Naveed et al.: Prediction of bioporeand matrix-dominated flow


Introduction
The importance of macropore flow for the partitioning of precipitation between runoff and infiltration, for plant water uptake and plant growth, for biogeochemical cycling rates, and for potential risks of ground water contamination, is widely recognized (Iversen et al., 2011;de Jonge et al., 2004;Fox et al., 2004;Moustafa, 2000).Thus, over the last decade, major research efforts have been devoted to improving the understanding of macropore flow and associated governing parameters, and to developing predictive macropore flow models (Jarvis, 2007).Macropore flow and transport refers to the localized and commonly very rapid movement of water and solutes through the soil profile.Macropores resulting from biological activity (root channels, worm holes, etc.), geological forces (subsurface erosion, shrinkage and swelling, etc.), and agricultural management (e.g., plowing) serve as the main channels for this rapid and long-distance flow and transport of water, air, and contaminants.Macropore flow is largely determined by soil structure and is generally a dominating process in loamy and clayey soils (Jarvis et al., 2009) where large inter-aggregate pores and biopores often act as pathways for rapid flow and transport.The transition from matrix to macropore flow (equilibrium to non-equilibrium) depends on the pore size distribution and pore continuity, and the degree of soil saturation (Bouma, 1981).Macropore flow often occurs in pores with equivalent effective cylindrical diameters larger than 0.3 mm, which indicates that the matric potential needs to be close to 0 and the water content close to saturation for these pores to be activated (Jarvis, 2007).
Soil and crop management practices strongly modify soil structure and thus the extent of macropore flow and transport.Wang et al. (2013) and Gonzalez-Sosa et al. (2010) studied the impact of land use on the hydraulic properties of the topsoil of the Loess Plateau of China and for a suburban catchment in France, respectively.Both studies have reported greater saturated hydraulic conductivities for forested land, intermediate for permanent pasture, and lower for farmland soils.This is primarily due to the abundance of biota and less disturbance in forests and permanent pastures when compared to cultivated lands (Naveed et al., 2014a;Norgaard et al., 2013;Pérèsa et al., 2012).Application of animal manure and fertilizers can also influence macropore flow, first by altering soil structure and second by promoting the density of earthworms, particularly deep penetrating anecic worms (Naveed et al., 2014b).Climatic conditions (seasonal temperature and precipitation variations) might also affect soil structure and macropore flow through interactions with physical processes such as cyclic freezing/thawing and wetting/drying (Hu et al., 2012).Due to the complex interactions and the significant number of influencing factors, a large spatial variability of saturated hydraulic conductivity has been reported by several authors (Wang et al., 2013;Raczkowski et al., 2012;Iversen et al., 2011).Therefore, the predictive capabilities of empirical models/pedotransfer functions for saturated hydraulic conductivity are limited because they ignore the effects of key site factors and underestimate the significance of soil structure (Vereecken et al., 2010).Recently, pedotransfer functions for saturated hydraulic conductivity that account for soil structure have been developed, but they are rarely applied due to the complexity of input parameters and the still significant prediction inaccuracies (Jarvis et al., 2013;Iversen et al., 2011;Lilly et al., 2008).
Along with the prediction of macropore water flow (i.e., saturated hydraulic conductivity), prediction of macropore airflow (i.e., air permeability and diffusivity) is also of essence.Air permeability is a key parameter for the design of soil vapor extraction remediation methods.Air diffusivity is of importance because the availability of oxygen to plant roots via diffusion is a basic factor for plant productivity.Various empirical models have been proposed in the past for the prediction of air permeability (Chamindu Deepagoda et al., 2011;Kawamoto et al., 2006) and air diffusivity (Chamindu Deepagoda et al., 2011;Moldrup et al., 2000).However, none of the above studies have evaluated their applicability after distinguishing between biopore-and matrixdominated flow domains.
Recent developments in soil imaging techniques not only allow visual observations, but also quantification of pore network complexity.Application of X-ray computed tomography (CT) provides emerging alternative means for estimating subsurface macropore flow and transport (Wildenschild and Sheppard, 2013).Over the last decade, numerous studies related to the characterization of macropore structure (i.e., macroporosity, macropore size distribution, volume, surface area, and tortuosity) have been conducted with X-ray CT for different land use and management systems (Katuwal et al., 2015;Larsbo et al., 2014;Hu et al., 2014;Naveed et al., 2013;Vogel et al., 2010;Luo et al., 2010).However, to date there have been only very few published studies on quantitatively relating macropore network characteristics to the observations of macropore flow.Katuwal et al. (2015) found that CT-derived macroporosity for the limiting section of a soil column was strongly correlated with air permeability and 5 % tracer arrival time.Larsbo et al. (2014) reported significant correlations between X-ray CT-derived macropore network characteristics and flow and transport parameters.Paradelo et al. (2013) found that CT-derived macroporosity was strongly correlated with saturated hydraulic conductivity, solute dispersivity, and contaminant breakthrough.Luo et al. (2010) reported that macroporosity, path number, hydraulic radius, and macropore angle were the most useful X-ray CT-derived parameters for predicting macropore flow and transport under saturated conditions.
In this study we first evaluate the predictive performance of existing pedotransfer functions/models for saturated hydraulic conductivity, air permeability, and gas diffusivity.
While it has been previously demonstrated that water flow in macropores cannot be accurately predicted with empirical models from basic soil properties (Weynants et al., 2009;Vereecken et al., 2010), there is only little published work related to gas diffusivity.Furthermore, existing pedotransfer functions/empirical models do not discern between matrixand biopore-dominated flow domains, which is of significance for understanding and accurate prediction of preferential flow as demonstrated in the results section.In the second part of this study we derive novel macropore network characteristics from X-ray CT observations for the prediction of saturated hydraulic conductivity, air permeability, and gas diffusivity, which demonstrated their utility for improving accuracy of gas and water flow predictions.The simplest form of the Kozeny-Carman equation proposed by Ahuja et al. (1984) is parameterized with novel CT-derived parameters such as effective percolating macroporosity for biopore-dominated flow and total macroporosity for matrixdominated flow, and improvement of prediction accuracy is discussed.

Study site and soil sampling
The 1.69 ha study site located in Silstrup in northwestern Denmark (56 • 55 56 N, 8 • 38 44 E) is composed of glacial till, a dominant geological formation covering about 43 % of all farmland in Denmark (Geological Survey of Denmark and Greenland, 1999).The top meter of the soil is highly fractured and bioturbated, containing 100 to 1000 biopores per m 2 .The field was not tilled for about 3 years prior to soil sampling.It has been plowed in December 2008 to 23 cm depth and harrowed twice to 5 cm depth in March 2009.Since then the soil was only disturbed when slurry was injected in 10 cm depth in April 2009 and in 5 cm depth in September 2009.A thorough overview of management practices at the study site between 2006 and 2010 is provided in Norgaard et al. (2013).
Sixty-five undisturbed cylindrical soil cores (6 cm inner diameter and 3.5 cm height) were extracted from the topsoil (5 to 8.5 cm depth) in summer 2012.At the time of sampling the field was cultivated with red fescue (Festuca rubra L.).The soil columns were sampled on a 15 m × 15 m grid with an additional five sampling locations between grid points (Fig. 2).All soil columns were extracted by pushing a customized core sampler with aluminum sampling cylinders into the soil and removing the surrounding material step by step.Extracted soil columns were immediately covered with tight plastic lids, placed in plastic bags, and carefully transported to the laboratory to avoid smearing and compaction effects.The soil columns were stored in an environmentally controlled room at 2 • C until the start of the measurements.In addition, bulk soil samples were collected from each point at the same soil depth for texture and organic carbon analysis.

X-ray computed tomography scanning and analysis
An industrial X-ray CT scanner (X-Tek HMX225) at the Helmholtz Center for Environmental Research in Halle in Germany was used to scan the intact soil columns at a voltage of 180 kV and a current of 400 µA.A copper filter was placed between the X-ray source and the soil columns to alleviate beam hardening.The shadow projections (radiographs) were reconstructed with a Feldkamp cone-beam algorithm (Feldkamp et al., 1984) to obtain 16-bit grayscale 3-D data comprised of 500 × 500 × 300 voxels at a resolution of 129 µm (Fig. 1a).For subsequent analysis, the 3-D grayscale volumes were cropped to remove the container wall and disturbed regions on the top and bottom of the sample, numerically corrected for intensity differences caused by beam hardening and other scanning artifacts with a sequential algorithm developed by Iassonov and Tuller (2010), and a 3-D median filter (Jassogne et al., 2007) with a radius of 6 voxels was applied to the grayscale volumes to remove noise (Fig. 1b).Though median filtering is computationally more demanding than conventional smoothing filters, it is less sensitive to outlier values, and thus preserves edges.A locally adaptive Bayesian Markov random field (MRF) algorithm (Iassonov et al., 2009;Kulkarni et al., 2012) that was seeded with adaptive K-means clustering (Chen et al., 1998) was used to segment the intensity-corrected and filtered data to distinguish macropores from the soil matrix (Fig. 1c).The homogeneity parameter β in the MRF model was set to 2. For details of the applied MRF segmentation algorithm, see Kulkarni et al. (2012) and Tuller et al. (2013).
The segmented CT data for each soil column were further analyzed with the Image-J software package (Rasband, 1997(Rasband, -2016) ) to obtain macroporosity, percolating macroporosity, effective percolating macroporosity, macropore specific surface area, macropore hydraulic radius, macropore mean diameter, macropore fractal dimension, macropore global connectivity, and macropore local connectivity (see the flowchart depicted in Fig. 1).Three-dimensional pore visualization was conducted with the Image-J plugin 3-D viewer.Based on 3-D visual observations, soil columns containing percolating biopores (round shape formed either by roots or earthworms) were separated and labeled as bioporedominated flow columns; the rest were labeled as matrixdominated flow columns (Fig. 1d).The number of pore voxels was determined from the segmented data, and macroporosity (MP) was then calculated as the ratio of the number of pore voxels to the number of total sample voxels (Fig. 1d).The percolating macroporosity (PMP) was calculated based on only the pores that were connected from sample top to bottom by removing all isolated pores (Fig. 1e).All isolated pores were removed with Image-J plugin "Find Connected Regions".Effective percolating macroporosity (EPMP) was defined and calculated as the ratio of the minimum cross-sectional area of percolating macropores (while moving voxel layer by voxel layer from the top to the bottom of the core) and the cross-sectional area of the soil column (Fig. 1f).Macropore specific surface area (MPSSA) was calculated as the ratio of the surface area of macropores and the volume of the soil column (Fig. 1g).This was accomplished with Image-J plugin "Analyze Particles".Macropore hydraulic radius (MPHR) was defined as the ratio of macropore volume and macropore surface area (Fig. 1h) applying Image-J plugin "Analyze Particles".The macropore mean diameter (MPMD) was estimated with a local 3-D thickness algorithm proposed by Dougherty and Kunzelmann (2007) and embedded in Image-J plugin "Bone-J".This algorithm defines the pore diameter as the diameter of the largest sphere that fits within the pore.The histogram of the thickness map was used for estimating macropore size distribution and macropore mean diameter (Fig. 1i).Macropore fractal dimension (MPFD) was calculated as a measure of the heterogeneity of the spatial distribution of macroporosity with Image-J plugin "Bone-J" (Fig. 1j).Macropore global connectivity (MPGC) was defined and calculated as the ratio of the percolating macroporosity and the total macroporosity of the soil column (Fig. 1k).The macropore local connectivity (MPLC) was estimated with Image-J plugin "Bone-J" (Fig. 1l).MPLC equals 1 if all pores are connected in one percolating cluster and 0 if porosity is fragmented into many clusters of similar size.X-ray CT-derived pore network characteristics for all scanned and analyzed core samples are provided in Supplement Table S1.

Soil physical measurements
Soil texture was determined from disturbed soil samples with a combination of wet sieving and the hydrometer method, after passing the sample through a 2 mm sieve.Soil organic carbon was determined with a LECO carbon analyzer (St.Joseph, MI, USA) coupled with an infrared CO 2 detector.A multiplication factor of 1.72 was used to convert soil organic carbon to soil organic matter.The sand, silt, clay, and organic matter contents for the 65 investigated samples are listed in Table S2.
After X-ray CT scanning, air permeability and gas diffusivity at −30 and −100 cm matric potentials, and saturated hydraulic conductivity (K sat ), were measured on the same columns.The soil columns were placed in a sand box and saturated from the bottom with tap water.After saturation, tension was successively applied to establish matric potentials of −30 and −100 cm, respectively.Air permeability (K a ) was then measured with the steady-state method described in Iversen et al. (2001) at both −30 and −100 cm matric potentials.A pressure of 5 hPa was applied to assure laminar flow during the measurements.The K a was calculated based on the Darcy equation considering the pressure difference across the soil cores: where ) is the dynamic viscosity of air, a s (L 2 ) is the cross-sectional area, and L s (L) is the length of the column.Gas diffusivities (D P /D 0 ) at −30 and −100 cm matric potentials were measured with the onechamber method developed by Schjønning et al. (2013).After D P /D 0 measurements, the soil columns were resaturated, and the saturated hydraulic conductivity (K sat ) was measured with the constant head method (Klute and Dirksen, 1986).All measured flow parameters are provided in Table S3.Ahuja et al. (1984) developed a relationship (EPM, effective porosity model) between saturated hydraulic conductivity (K sat ) and effective porosity (φ e ) based on the generalized Kozeny-Carman equation:

Modeling
where K sat is saturated hydraulic conductivity, K a is air permeability, D P /D 0 is gas diffusivity, and A and B are empirical constants.Ahuja et al. (1984) defined φ e as the total porosity minus the soil volumetric water content at field capacity assumed at a matric potential of −33 kPa.Based on a simple calculation applying the capillary rise equation, this means that φ e is the porosity contributed by pores larger than about 9 µm in diameter.We first parameterized the original Ahuja et al. (1984) model with φ e equivalent to the air-filled porosity at −30 kPa.Then, X-ray CT-derived macroporosity (MP) was used for φ e for matrix-dominated flow, and X-ray CT-derived effective percolating macroporosity (EPMP) was applied for φ e for biopore-dominated flow.Note that because of the 129 µm resolution of the CT scans, the CT-derived parameters MP and EPMP represent significantly larger pores than originally suggested in Ahuja et al. (1984).This seems quite reasonable and interesting to test as macropore flow often occurs in pores with equivalent effective cylindrical diameters larger than 300 µm (Jarvis, 2007).Rawls et al. (1998) reported that several researchers found the slope A to vary between 1.59 and 3.98 and the intercept to vary between 440 and 34 000 cm d −1 .

Statistics
Data collected for soil textural properties and macropore flow parameters were first subjected to classical statistical analysis to obtain descriptive statistics, including minimum, maximum, mean, median, standard deviation, skewness, and coefficient of variation (CV).The degree of spatial variability of soil textural properties and macropore flow parameters was determined with ordinary Kriging.The ArcMap 10.1 software (Esri Inc., Redlands, CA, USA) was used to generate contour maps for each measured soil property.Spearman rank order correlation coefficients between macropore network characteristics and macropore flow parameters were calculated with the commercial SigmaPlot 11.0 software package (Systat Software, Inc., San Jose, CA, USA).Selected correlations were also graphically displayed and analyzed with linear, power, or exponential regression models.While the applicability of linear models was evaluated, power or exponential models yielded significantly better results in most cases.The models were only fitted if they were significant at p < 0.01.

Results and discussion
3.1 Spatial variability of soil texture, organic matter, and macropore flow parameters The soil of the study site was classified as sandy loam (USDA-NRCS Web Soil Survey, 2010) with clay contents ranging from 14 to 19 %, and organic matter contents varying from 2.9 to 3.8 %.Descriptive statistics for all soil textural properties are depicted in Table 1.Clay and sand contents were positively skewed, whereas silt and organic matter contents were negatively skewed.All soil textural properties were slightly variable across the field, with coefficients of variation (CV) below 10 %.It has been previously reported that the CV for soil textural properties generally depends on the extent of the study area.For example, Sharma et al. (2011) reported a CV for soil textural properties within the range of 20 to 30 % for a 40 ha agricultural field in New Mexico, while Wang et al. (2013) reported a CV within the range of 19 to 156 % across the Loess Plateau of China (620 × 10 3 km 2 ).Kriged maps indicated that soils with high clay contents (Fig. 2a) were on the northern side of the field, whereas soils with high organic matter contents occupied the southern side (Fig. 2d).Thus, clay and organic matter gradients run in opposite directions at the study site.Soils with high silt contents (Fig. 2b) were in the western part of the field, whereas soils with high sand contents were in the eastern part (Fig. 2c).Relevant information about the semivariograms for each interpolated map is provided in Table 2.
Descriptive statistics for saturated hydraulic conductivity (K sat ), air permeability (K a ), and gas diffusivity (D P /D 0 ) at −30 and −100 cm matric potentials are provided in Table 1.Large positive skewness and quite different mean and median values were observed for all five macropore flow parameters.The K sat , K a , and D P /D 0 at −30 and −100 cm matric potentials showed the largest variations across the study site with a CV ranging from 92 to 218 %.High CV values were observed due to the presence of biopores in some of the soil columns, while not in others.Renderings of the samples marked as I, II, III, and IV in Fig. 2 are depicted in Fig. 3. Samples I and II are matrix-flow-dominated and samples III and IV are biopore-flow-dominated.Irrespective of the extent of the study area, large variations in K sat were also reported in other studies (e.g., Wang et al., 2013;Sharma et al., 2011;and Iqbal et al., 2005).Kriged maps for K sat , K a , and D P /D 0 (Fig. 2e-g) look quite similar, with some areas randomly exhibiting a high level of macropore flow, while matrix flow dominated in other regions irrespective of soil texture and organic matter content.

Predictive performance of empirical models
For many hydrological applications, saturated hydraulic conductivity (K sat ) is estimated from more readily available proxy variables such as texture and bulk density.Various empirical models/pedotransfer functions (e.g., Iversen et al., 2011;Jarvis et al., 2009;Schaap et al., 2001;Wösten et al., 1999;Revil and Cathles, 1999) have been previously proposed for predicting saturated hydraulic conductivity.We have observed poor predictive performance of empirical K sat models such as proposed by Revil and Cathles (1999) and Schaap et al. (2001) (Fig. 4) and for models proposed by Wösten et al. (1999), Vereecken et al. (1989), and Cosby et al. (1984) (not shown).While the measured saturated hydraulic conductivities span over 5 orders of magnitude due to the presence of a wide range of macropores and biopores in the core samples, model predictions yielded a very narrow K sat range (Fig. 4).The primary reason for the failure of existing empirical models/pedotransfer functions is that they only consider soil texture and bulk density, and thus are not able to realistically capture macropore flow, particularly for highly structured and bioturbated soils.In general, empirical models overpredicted K sat in the case of matrix flow (empty symbols), while they underpredicted K sat for soil columns with biopore flow (filled symbols).Because results were obtained for samples of limited size from the A horizon, it  should be noted that for larger scales the structural characteristics and associated flow parameters, especially the parameters related to pore connectivity, might change.

Correlations between macropore flow parameters and macropore network characteristics
The CT-derived macroporosity and the physically measured air-filled porosity at −30 cm matric potential are in good agreement as shown in Fig. 6.At −30 cm matric potential, all pores with diameters larger than 100 µm should have drained according to the capillary-rise equation.This indicates that the physically measured air-filled porosity at −30 cm matric potential (pore diameter > 100 µm) should be greater than the X-ray CT-derived macroporosity (resolution = 129 µm).However, this is only true when assuming a parallel bundle of capillary tubes, which is not a realistic assumption for natural soils.Due to the ink-bottle effect a considerable number of pores with diameters > 100 µm are expected to be water filled after drainage at a matric potential of −30 cm.Hence, no perfect match between the CT-measured morphological pore size and the hydraulic pore size estimated with the capillary-rise equation should be expected (Vogel, 2000).The observed agreement between the two measures is reasonable and confirms the applicability of the applied image segmentation method (Fig. 6).However, it must be noted that different image segmentation methods can result in quite different macroporosity values if the CT image quality is bad,  (c, d) WLR-Marshall model (Moldrup et al., 2000).Filled symbols represent samples with biopore-dominated flow and empty symbols represent samples with matrix-dominated flow.Visualizations of samples marked as I, II, III, and IV are depicted in Fig. 3. i.e., if there is a lot of noise and partial volume effects as shown in Naveed (2014).Spearman rank order correlation analysis for macropore flow parameters and macropore network characteristics was performed for all soil columns (Fig. 7a), for soil columns with biopore(s) connected from the top to the bottom (Fig. 7b), and for soil columns with inter-aggregate macropores or disconnected biopores (Fig. 7c).Many of the CT-derived macropore network characteristics were strongly correlated (Fig. 7).This is because large macroporosity is associated with large macropore surface area and better connectivity of macropores.This is in agreement with other recent studies (e.g., Katuwal et al., 2015, andLarsbo et al., 2014).The macropore mean diameter and hydraulic radius were however poorly correlated with other macropore network characteristics.Significant Spearman rank order cor- MP is macroporosity, PMP is percolating macroporosity, EPMP is effective percolating macroporosity, MPSSA is macropore specific surface area, MPHR is macropore hydraulic radius, MPMD is macropore mean diameter, MPFD is macropore fractal dimension, MPGC is macropore global connectivity, MPLC is macropore local connectivity, K sat is saturated hydraulic conductivity (cm h −1 ), K a − 30 is air permeability (µm 2 ) at −30 cm matric potential, K a − 100 is air permeability (µm 2 ) at −100 cm matric potential, D P /D 0 − 30 is gas diffusivity at −30 cm matric potential, and D P /D 0 − 100 is gas diffusivity at −100 cm matric potential.Strong correlation (r > 0.70), moderate correlation (r = 0.5-0.7),and weak correlation (r < 0.5).relations were also observed between macropore flow parameters and most of the CT-derived macropore network characteristics (Fig. 7).X-ray CT-derived macroporosity was strongly correlated with macropore flow parameters for all three categories of soil samples (Fig. 7a, b, and c).Very strong correlations were observed between effective percolating macroporosity (EPMP) and macropore flow parameters for the soil columns with biopores connected from Table 3. Empirical constants for the Ahuja (1984) model with air-filled porosity at −30 kPa, X-ray CT-derived effective percolating macroporosity (EPMP), and total macroporosity (MP) as effective porosity φ e , respectively.

Variable
A B φ e : air-filled porosity at −30 kPa in the Ahuja (1984) model Saturated hydraulic conductivity, K sat (cm h −1 ) 5000 3.2 Air permeability at −30 cm, K a − 30 (µm 2 ) 5000 3.4 Air permeability at −100 cm, K a − 100 (µm 2 ) 5000 3.2 Gas diffusivity at −30 cm, D P /D 0 − 30 0.27 2.3 Gas diffusivity at −100 cm, D P /D 0 − 100 0.27 2.0 φ e : effective percolating macroporosity (EPMP) in the Ahuja (1984) model for biopore-dominated flow Saturated hydraulic conductivity, K sat (cm h −1 ) 5000 1.15 Air permeability at −30 cm, K a − 30 (µm 2 ) 5000 1.5 Air permeability at −100 cm, K a − 100 (µm 2 ) 5000 1.4 Gas diffusivity at −30 cm, D P /D 0 − 30 0.27 1.12 Gas diffusivity at −100 cm, D P /D 0 − 100 0.27 0.98 φ e : total macroporosity (MP) in the Ahuja (1984) model for matrix-dominated flow Saturated hydraulic conductivity, K sat (cm h −1 ) 5000 3.2 Air permeability at −30 cm, K a − 30 (µm 2 ) 5000 3.0 Air permeability at −100 cm, K a − 100 (µm 2 ) 5000 2.7 Gas diffusivity at −30 cm, D P /D 0 − 30 0.27 1.90 Gas diffusivity at −100 cm, D P /D 0 − 100 0.27 1.55 the top to the bottom (Fig. 7b).Macropore hydraulic radius and macropore mean diameter were significantly correlated with macropore flow parameters for the soil columns associated with biopore-dominated flow (Fig. 7b), but poorly correlated for soil columns associated with matrix-dominated flow (Fig. 7c).These findings are in agreement with Elliot et al. (2010) and Quinton et al. (2008).Both macropore global and local connectivity were poorly correlated with macropore flow parameters for the soil columns associated with biopore-dominated flow (Fig. 7b), but significantly correlated for the soil columns associated with matrix-dominated flow (Fig. 7c).This makes sense as biopore flow is mainly governed by the largest biopore present in the soil column, whereas matrix flow is mainly controlled by the pore size distribution and connectivity of pores.Selected correlations were graphically displayed and analyzed with linear, power, and exponential regression models.The latter were superior to linear models in most cases, as shown in Fig. 8.The saturated hydraulic conductivity (K sat ) was plotted as a function of CT-derived macroporosity (Fig. 8a).Two distinct branches were observed for lower macroporosity values, which approach a single branch with increasing CT-derived macroporosity.The upper branch with greater conductivities comprises core samples with one or more biopores connected from top to bottom that mainly govern fluid flow (filled symbols).Samples III and IV marked in Fig. 8a and shown in Fig. 3 are members of this branch.The lower branch consists of core samples with fluid mainly flowing through inter-aggregate and textural pores (empty symbols).Samples I and II marked in Fig. 8a and shown in Fig. 3 are members of this branch.Distinct significant power regressions were observed between K sat and macroporosity for these two categories of the soil columns (Fig. 8a).This suggests that biopore-dominated and matrix-dominated flow columns should be discerned as an initial step prior to studying the relationships between macropore flow and CT-derived macroporosity.Both Paradelo et al. (2013) and Luo et al. (2010) found similar relationships between saturated hydraulic conductivity and CTderived macroporosity, with R 2 ranging from 0.50 to 0.60.A stronger power regression was observed when K sat was plotted as a function of the effective percolating macroporosity (R 2 increased from 0.43 to 0.76), for the soil columns associated with biopore-dominated flow (Fig. 8b, filled symbols), but this is not the case for the soil columns with matrix-dominated flow (Fig. 8b, empty symbols).Significant power regressions were observed between K sat and macropore mean diameter (Fig. 8c).No significant regression was observed between K sat and macropore local connectivity as shown in Fig. 8d.No significant regression was observed between K sat and macropore local connectivity for the soil samples associated with biopore-dominated flow (Fig. 8d, filled symbols).This suggests that connectivity of the macropores is not the controlling factor in case of K sat .
Air permeability at −30 cm matric potential, K a (−30), was plotted as a function of macroporosity as shown in Figure 8. Saturated hydraulic conductivity (K sat ), air permeability at −30 cm matric potential (K a − 30), air permeability at −100 cm matric potential (K a − 100), gas diffusivity at −30 cm matric potential (D P /D 0 − 30), and gas diffusivity at −100 cm matric potential (D P /D 0 − 100) were plotted as a function of selected CT-derived macropore network characteristics; filled symbols represent samples with biopore-dominated flow and empty symbols represent samples with matrix-dominated flow.Fitting linear regression models has been attempted; a power/exponential model was always superior where a significant correlation was present.Two separate regressions were fitted for samples with biopore flow and matrix flow if they were significantly different.Plots (g), (k), (l), and (p) only show one curve because the other was not significant, while plots (q), (r), and (t) have only one model because the two models did not differ significantly from each other.Fig. 8e.Distinct significant power regressions were observed for the two categories of soil columns, i.e., columns with biopore-dominated flow and with matrix-dominated flow (Fig. 8e).Similar to K sat , the power regression was significantly improved (R 2 increased from 0.49 to 0.80) when K a (−30) was plotted as a function of effective percolating macroporosity instead of total macroporosity for the soil columns associated with biopore-dominated flow (Fig. 8f, filled symbols).A significant power regression was observed between K a (−30) and macropore mean diameter for the soil columns with biopore-dominated flow, while no significant regression was observed between K a (−30) and macropore mean diameter for the soil columns with matrix-dominated flow (Fig. 8g).No significant regression was observed between K a (−30) and macropore local connectivity (Fig. 8h).Similar power regressions were also observed for K a (−100) as a function of macroporosity, effective percolating macroporosity, and macropore mean diameter as shown in Figs. 8i, j, Hydrol. Earth Syst. Sci., 20, 4017-4030, 2016 www.hydrol-earth-syst-sci.net/20/4017/2016/ and k, respectively.The K a (−100) showed significant exponential regression as a function of macropore local connectivity for matrix dominated-flow columns (Fig. 8l). Figure 8m and n showed significant power regressions when gas diffusivity at −30 cm matric potential, D P /D 0 (−30), was plotted against macroporosity and effective percolating macroporosity, respectively.Distinct significant power regressions observed for soil columns associated with biopore-dominated flow and matrix-dominated flow reflect the fact that preferential diffusive flow occurred at −30 cm matric potential.However, at −100 cm matric potential, a single regression significantly described both types of data associated with biopore flow and matrix flow as shown in Fig. 8q and r.This indicates that no preferential diffusive flow occurred at and below −100 cm matric potentials.Both D P /D 0 (−30) and D P /D 0 (−100) showed insignificant regressions when plotted as a function of macropore mean diameter for both categories of soil samples (Fig. 8o and  8s).Significant exponential regressions were observed when D P /D 0 (−30) and D P /D 0 (−100) were plotted as a function of macropore local connectivity for both sets of soil columns associated with matrix flow and biopore flow (Fig. 8p and  t).This is expected as D P /D 0 is a concentration-driven gas transport parameter mainly controlled by total air-filled pore space and its connectivity, and not by the pore size (Moldrup et al., 2000).

Modeling saturated hydraulic conductivity, air permeability, and diffusivity
Saturated hydraulic conductivity, air permeability at −30 and −100 cm matric potentials, and gas diffusivity at −30 and −100 cm matric potentials were modeled with the simplified form of the Kozeny-Carman equation presented in Ahuja et al. (1984).First, we have tested the predictive performance of the original Ahuja et al. (1984) model with air-filled porosity at −30 kPa as the effective porosity (Fig. 9, red empty symbols).Then, we have modified the original equation with novel CT-derived input parameters.The effective porosity in the original model was replaced with the CT-derived total macroporosity (MP) in case of matrix-dominated flow (Fig. 9, black empty symbols) and with the effective percolating macroporosity (EPMP) in case of biopore-dominated flow (Fig. 9, black filled symbols).The empirical fitting parameters (A and B) for saturated hydraulic conductivity, air permeability at −30 and −100 cm matric potentials, and gas diffusivity at −30 and −100 cm matric potentials are provided in Table 3.The 1 : 1 plots of measured and predicted saturated hydraulic conductivity, air permeability, and gas diffusivity are shown in Fig. 9. From Fig. 9 it is obvious that predictions with the Ahuja et al. (1984) model with novel input data from X-ray CT analysis are very reasonable and yielded better results than the conventionally parameterized Ahuja et al. (1984) model.This indicates that Xray CT-derived macropore characteristics (MP and EPMP) at 129 µm resolution are quite useful for predicting macropore flow.However, distinguishing between biopore-and matrixdominated flows is a prerequisite.The predictive capability of the proposed modeling framework requires further independent validation for different soil types to confirm the values/ranges for empirical constants A and B for saturated hydraulic conductivity, air permeability, and gas diffusivity.

Conclusions and perspective
While soil textural properties exhibited small spatial variability across the study site with a CV < 10 %, the macropore flow parameter saturated hydraulic conductivity, air permeability, and gas diffusivity showed large spatial variability across the field with a CV > 100 %.The predictive performance of existing empirical models/pedotransfer functions for saturated hydraulic conductivity and air permeability at −30 and −100 cm matric potentials was unsatisfactory.For saturated hydraulic conductivity, existing empirical models overpredicted for cases with matrix-dominated flow and underpredicted for cases with biopore-dominated flow.With regard to air permeability, em-pirical models predicted matrix-dominated flow reasonably well, whereas significant underpredictions were observed for cases with biopore-dominated flow.The tested empirical model for the prediction of gas diffusivity performed well at −100 cm matric potential, while it failed at −30 cm matric potential, particularly for the soil columns that contained biopores that were connected from the sample top to the sample bottom (i.e., biopore-flow-dominated samples).
Significant Spearman rank correlations were observed between CT-derived macropore network characteristics and macropore flow parameters.These correlations were further improved when the soil columns were separated into matrixdominated and biopore-dominated flow columns.The predictive performance of the Ahuja et al. (1984) model with novel input parameters, namely X-ray CT-derived effective percolating macroporosity (EPMP) for biopore-dominated flow and total macroporosity (MP) for matrix-dominated flow, was significantly improved.However, further studies for different soil textures are required to confirm the values/ranges of the empirical Ahuja et al. (1984) A and B model parameters for accurate predictions of saturated hydraulic conductivity, air permeability, and gas diffusivity.
The rapid development of advanced CT-image segmentation and analysis tools in conjunction with computational fluid dynamics provides promising future means to simulate the dynamics of flow and transport directly with CT-derived macropore networks as boundaries.One method particularly suitable for simulating macropore flow and transport based on X-ray CT data is the lattice Boltzmann method (LBM).Most of the studies to date that applied the LBM for simulating flow and transport based on CT data were for granular porous media (glass beads/sand) and fractured rocks, and not for natural field soil samples.The strong correlations between macropore flow parameters and X-ray CTderived macropore network characteristics observed in this study suggest that lattice Boltzmann flow and transport simulations based on X-ray CT images are a promising avenue for future research.

Figure 1 .
Figure 1.Flowchart illustrating all performed CT-data enhancement, segmentation, and analysis steps.

Figure 3 .
Figure 3. Three-dimensional visualizations of sample soil columns and associated measured macropore flow parameters.K sat is saturated hydraulic conductivity, and K a −100 and D P /D 0 −100 are air permeability and gas diffusivity at −100 cm matric potential, respectively.

Figure 4 .
Figure 4. Predictive performance of empirical saturated hydraulic conductivity (K sat ) models.Filled symbols represent samples with biopore-dominated flow and empty symbols represent samples with matrix-dominated flow.Visualizations of samples marked as I, II, III, and IV are depicted in Fig. 3.

Figure 5 .
Figure 5. Predictive performance of empirical models for air permeability (K a ) and gas diffusivity (D P /D 0 ) at −30 and −100 cm matric potentials.(a, b) Chamindu Deepagoda et al. (2011) model;(c, d) WLR-Marshall model(Moldrup et al., 2000).Filled symbols represent samples with biopore-dominated flow and empty symbols represent samples with matrix-dominated flow.Visualizations of samples marked as I, II, III, and IV are depicted in Fig.3.

Figure 6 .
Figure 6.CT-derived macroporosity plotted as a function of physically measured air-filled porosity at −30 cm matric potential.

Table 1 .
Descriptive statistics for selected soil physical properties (no. of samples = 65).

Table 2 .
Partial sill, nugget, range, Kriging interpolation model, and root mean square error (RMSE) for semivariograms for each interpolated map.All interpolations were performed with ESRI ArcMap 10.1.