Saturated hydraulic conductivity model computed from bimodal water retention curves for a range of New Zealand soils

Descriptions of soil hydraulic properties, such as the soil moisture retention curve, θ(h), and saturated hydraulic conductivities, Ks, are a prerequisite for hydrological models. Since the measurement of Ks is expensive, it is frequently derived from statistical pedotransfer functions (PTFs). Because it is usually more difficult to describe Ks than θ(h) from pedotransfer functions, Pollacco et al. (2013) developed a physical unimodal model to compute Ks solely from hydraulic parameters derived from the Kosugi θ(h). This unimodal Ks model, which is based on a unimodal Kosugi soil pore-size distribution, was developed by combining the approach of Hagen–Poiseuille with Darcy’s law and by introducing three tortuosity parameters. We report here on (1) the suitability of the Pollacco unimodal Ks model to predict Ks for a range of New Zealand soils from the New Zealand soil database (S-map) and (2) further adaptations to this model to adapt it to dual-porosity structured soils by computing the soil water flux through a continuous function of an improved bimodal pore-size distribution. The improved bimodal Ks model was tested with a New Zealand data set derived from historical measurements of Ks and θ(h) for a range of soils derived from sandstone and siltstone. The Ks data were collected using a small core size of 10 cm diameter, causing large uncertainty in replicate measurements. Predictions of Ks were further improved by distinguishing topsoils from subsoil. Nevertheless, as expected, stratifying the data with soil texture only slightly improved the predictions of the physical Ks models because the Ks model is based on pore-size distribution and the calibrated parameters were obtained within the physically feasible range. The improvements made to the unimodal Ks model by using the new bimodalKs model are modest when compared to the unimodal model, which is explained by the poor accuracy of measured total porosity. Nevertheless, the new bimodal model provides an acceptable fit to the observed data. The study highlights the importance of improving Ks measurements with larger cores.

Abstract.Descriptions of soil hydraulic properties, such as the soil moisture retention curve, θ (h), and saturated hydraulic conductivities, K s , are a prerequisite for hydrological models.Since the measurement of K s is expensive, it is frequently derived from statistical pedotransfer functions (PTFs).Because it is usually more difficult to describe K s than θ (h) from pedotransfer functions, Pollacco et al. (2013) developed a physical unimodal model to compute K s solely from hydraulic parameters derived from the Kosugi θ (h).This unimodal K s model, which is based on a unimodal Kosugi soil pore-size distribution, was developed by combining the approach of Hagen-Poiseuille with Darcy's law and by introducing three tortuosity parameters.We report here on (1) the suitability of the Pollacco unimodal K s model to predict K s for a range of New Zealand soils from the New Zealand soil database (S-map) and (2) further adaptations to this model to adapt it to dual-porosity structured soils by computing the soil water flux through a continuous function of an improved bimodal pore-size distribution.The improved bimodal K s model was tested with a New Zealand data set derived from historical measurements of K s and θ(h) for a range of soils derived from sandstone and siltstone.The K s data were collected using a small core size of 10 cm diameter, causing large uncertainty in replicate measurements.Predictions of K s were further improved by distinguishing topsoils from subsoil.Nevertheless, as expected, stratifying the data with soil texture only slightly improved the predictions of the physical K s models because the K s model is based on pore-size distribution and the calibrated parameters were obtained within the physically feasible range.The improvements made to the unimodal K s model by using the new bi-modal K s model are modest when compared to the unimodal model, which is explained by the poor accuracy of measured total porosity.Nevertheless, the new bimodal model provides an acceptable fit to the observed data.The study highlights the importance of improving K s measurements with larger cores.

Introduction
Modelling of the water budget, irrigation, and nutrient and contaminant transport through the unsaturated zone requires accurate soil moisture retention, θ (h), and unsaturated hydraulic conductivity, K(θ ), curves.The considerable time and cost involved in measuring θ (h) and K(θ ) directly for a range of soils mean that the information for specific soils of interest is often not available (Webb, 2003).Therefore, these curves are generally retrieved from pedotransfer functions (PTFs), which are statistical relationships that generate lower-precision estimates of physical properties of interest based on many rapid and inexpensive measurements (e.g.Balland and Pollacco, 2008;Pollacco, 2008;Anderson and Bouma, 1973;Webb, 2003;Cichota et al., 2013).
The S-map database (Lilburne et al., 2012;Landcare Research, 2015) provides soil maps for the most intensively used land in New Zealand and is being gradually extended to give national coverage.S-map provides data for extensively used soil models, such as the soil nutrient model OVERSEER and the daily simulation model APSIM used by agricultural scientists.McNeill et al. (2012) used the New Zealand National Soils Database to derive PTFs to es-Published by Copernicus Publications on behalf of the European Geosciences Union.timate θ (h) at five tensions from morphological data of soils mapped in S-map.One of the current weaknesses of S-map is a lack of capacity to estimate K(θ ).Building on the work of Griffiths et al. (1999), Webb (2003) showed that morphologic descriptors for New Zealand soils can be used to predict K s .However, the predictions of K s were found to be too coarse for application to the wide range of soils within Smap.Therefore, Cichota et al. (2013) tested published statistical PTFs developed in Europe and the USA to predict θ (h) and K(θ ) for a range of New Zealand soils.They combined the best two or three PTFs to construct ensemble PTFs.They considered the ensemble PTF for θ (h) to be a reasonable fit, but the ensemble PTF for estimating K s exhibited large scatter and was not as reliable.The poor performance when estimating K s was possibly due to the absence of any measurements of pore-size distribution in their physical predictors (Watt and Griffiths, 1988;McKenzie and Jacquier, 1997;Chapuis, 2004;Mbonimpa et al., 2002) and also to the large uncertainties in the measurements from small cores (McKenzie and Cresswell, 2002;Anderson and Bouma, 1973).Consequently, there is an urgent need in New Zealand to develop a physically based K s model which is based on pore-size distribution.
Since PTFs developed to characterize θ (h) are more reliable than PTFs to characterize K(θ ) (e.g.Balland and Pollacco, 2008;Cichota et al., 2013), Pollacco et al. (2013) developed a new physical model that predicts unimodal K s solely from hydraulic parameters derived from the Kosugi (1996) θ (h).The K s model is derived by combining the Hagen-Poiseuille and Darcy laws (Sutera and Skalak, 1993) and by incorporating three semi-empirical tortuosity parameters.The model is based on the soil pore-size distribution and has been successfully validated using the European HYPRES (Wösten et al., 1998(Wösten et al., , 1999;;Lilly et al., 2008) and the UNSODA databases (Leij et al., 1999;Schaap and van Genuchten, 2006) but has not yet been applied to New Zealand soils.Most New Zealand soils are considered to be structured, with two-stage drainage (Carrick et al., 2010;McLeod et al., 2008) and bimodal pore-size distribution (e.g.Durner, 1994).Romano and Nasta (2016) showed by using the HYDRUS-1D package that large errors arise in the computation of the water fluxes if unimodal θ (h) and K(θ ) are used in structured soils.We therefore propose to improve the unimodal Pollacco et al. (2013) K s model so that it can predict K s for structured soils with bimodal porosity.
Measured K s values are widely recognized as one of the most variable soil attributes (McKenzie and Cresswell, 2002;Carrick, 2009).This is also recognized for New Zealand soils, both due to the high variability over short distances in soil parent material, age, depth, and texture, as well as strong macropore development with preferential macropore flow recognized as the norm rather than the exception in New Zealand soils (Webb et al., 2000;Carrick, 2009;McLeod et al., 2008).The measurement variability is also expected to increase as the sampling diameter decreases because small cores provide an unrealistic representation of the abundance and connectivity of macropores (McKenzie and Cresswell, 2002;Anderson and Bouma, 1973).McKenzie and Cresswell (2002) suggest that the standard Australian laboratory measurements should use cores with a minimum diameter of 25 cm and length of 20 cm.In New Zealand, K s has been obtained by using small cores, commonly with 10 cm diameter and 7.5 cm length.This has contributed to very high variability in measured K s (Webb et al., 2000).
The objectives of this research were to test the suitability of the unimodal Pollacco et al. (2013) K s model to predict K s from New Zealand soils, develop a K s bimodal model that makes predictions in structured soils solely from hydraulic parameters derived from the Kosugi θ (h), derive the uncertainties of the predictions of the K s bimodal model, and provide recommendations on the critical data sets that are required to improve the S-map database in New Zealand.

Kosugi unimodal water retention and unsaturated hydraulic conductivity curve
There are a number of closed-form unimodal expressions in the literature that compute the soil moisture retention curve θ (h) and the unsaturated hydraulic conductivity K(θ ) curves, such as the commonly used van Genuchten (1980) and Brooks and Corey (1964) curves.We selected the physically based Kosugi (1996) closed-form unimodal log-normal function expression of θ (h) and K(θ ) because its parameters are theoretically sound and relate to the soil pore-size distribution (Hayashi et al., 2009).Soils have a large variation in pore radius, r, which follows a log-normal probability density function.The unimodal Kosugi log-normal probability density function of pore radius (r) is often written in the following form: where θ r and θ s (cm 3 cm −3 ) are the residual and saturated water contents, r m (cm) is the median pore radius, and σ (−) denotes the standard deviation of ln(r).
Let S e denote the effective saturation, defining S e (r) = (θ − θ r ) / (θ r − θ s ) such that 0 ≤ S e ≤ 1. Integrating Eq. (1) from 0 to r yields the unimodal water retention curve as a function of r: where erfc is the complementary error function.The Young-Laplace capillary equation relates the soilpore radius, r, to the equivalent matric suction head, h (cm), at which the pore is filled or drained (i.e.r = Y / h, where Y = 0.149 cm 2 ).Kosugi's unimodal moisture retention curve θ uni (h) can be written in terms of S e : where h m (cm) is the median metric head.
The unimodal Kosugi unsaturated hydraulic conductivity function K(θ ) is written as where K s (cm day −1 ) is the saturated hydraulic conductivity.θ s is computed from the total porosity, φ, which is deduced from bulk density (ρ b ) and soil particle density (ρ p ) as follows: Due to air entrapment, θ s seldom reaches saturation of the total pore space φ (Carrick et al., 2011).Therefore, to take into account the fact that not all pores are connected, we perform the following correction of φ with α in the range [0.9, 1]: It is accepted that α = 0.95 (Rogowski, 1971;Pollacco et al., 2013;Haverkamp et al., 2005;Leij et al., 2005), but in this study the optimal α was found to be 0.98, since using a value of 0.95 resulted in several soil samples with θ 5 (θ measured at 5 kPa) greater than θ s , which is not physically plausible.This was due to the inaccuracy of measuring φ (discussed in Sect.4.1).The feasible range of the Kosugi hydraulic parameters is summarized in Table 1.The h m and σ feasible range is taken from Pollacco et al. (2013), who combined data from the HYPRES (Wösten et al., 1998(Wösten et al., , 1999;;Lilly et al., 2008) and UNSODA (Leij et al., 1999;Schaap and van Genuchten, 2006) databases.

Pollacco unimodal saturated hydraulic conductivity model
The saturated hydraulic conductivity model, K s_uni (Pollacco et al., 2013), computes K s from the Kosugi parameters θ s , θ r , σ , and h m (or r m ).K s_uni is based on the pore-size distribution (Eq. 1) and the tortuosity of the pores.K s_uni was derived by adopting the method of Childs and Collisgeorge (1950) and modelling the soil water flux through a continuous function of Kosugi (1996) pore-size distribution.This was performed by combining the Hagen-Poiseuille equation (Sutera and Skalak, 1993) with Darcy's law and introducing the connectivity and tortuosity parameters τ 1 , τ 2 of Fatt and Dykstra (1951) and τ 3 of Vervoort and Cattle (2003).K s_uni is computed as with for water at 20 • C, density of water ρ w = 0.998 g cm −3 , acceleration due to gravity g = 980.66cm s −2 , dynamic viscosity of water η = 0.0102 g cm −1 s −1 , and C is a constant equal to 1.03663 × 10 9 cm day −1 .
Integrating with S e instead of r avoids the complication of finding the minimum and maximum values of r.Isolating r of Eq. ( 2b) and replacing it in Eq. ( 7) gives dS e (8a) or and r m = Y /h m (Young-Laplace capillary equation) where τ 1 , τ 2 , τ 3 are tortuosity parameters [0-1).If tortuosity were not included (τ 1 , τ 2 , τ 3 = 0), the poresize distribution model would mimic the permeability of a bundle of straight capillary tubes.Vervoort and Cattle (2003) state "In reality soils are much more complex, with twisted and crooked pores, dead-ending or connecting to other pores.This means that there is a need to scale the permeability from the capillary tube model to include increased path length due to crookedness of the path (tortuosity) or lack of connection between points in the soil (connectivity)".Soils that are This takes into account the increased path length due to crookedness of the path.When τ 1 = 0, the flow path is perfectly straight down.When τ 1 increases, the flow path is no longer straight but meanders.
τ 2 This theoretically represents the shape of a microscopic capillary tube.The τ 2 parameter is used to estimate restrictions in flow rate due to variations in pore diameter and pore shape.When τ 2 = 0, the shape of the capillary tube is perfectly cylindrical.When τ 2 increases, the tube becomes less perfectly cylindrical, which causes lower connectivity.
τ 3 High-porosity soils tend to have large effective pores, θ s − θ r , which tend to be more connected than soils with smaller effective pores, which have more dead ends.When τ 3 = 0, the connectivity is the same between high-and low-porosity soils.When τ 3 increases, the connectivity of the soil increases (Vervoort and Cattle, 2003;Pollacco et al., 2013).Pollacco et al. (2013) found τ 3 to be the least sensitive parameter.

Romano bimodal water retention curve
New Zealand soils are predominantly well structured, with two-stage drainage (Carrick et al., 2010;McLeod et al., 2008) and therefore have a bimodal pore-size distribution (e.g.Durner, 1994).As K s_uni is based on a unimodal curve, θ uni (h), the proposed bimodal model, K s_bim , should be based on a bimodal θ bim (h) curve.Borgesen et al. (2006) showed that structured soils have both matrix (inter-aggregate) pore spaces and macropore (intra-aggregate) pore spaces.Thus, when the pores are initially saturated, such as (r > R mac ) or (h < H mac ), the flow is considered macropore flow, and when the soil is desaturated, such as (r < R mac ) or (h > H mac ), the flow is considered matrix flow, as shown in Fig. 1.R mac is the theoretical pore-size r that delimits macropore and matrix flow, and H mac is the theoretical pressure that delimits macropore and matrix flow.To model bimodal pore-size distribution, Durner (1994) superposes two unimodal pore-size distributions by using an empirical weighting factor, W , which partitions the volumetric percentage of macropore and matrix pores.Recently, Romano et al. (2011) proposed the following Kosugi bimodal θ bim_rom (h) distribution: where θ s , h m_mac , and σ _mac are, respectively, the saturated water content, the median pore radius, and the standard deviation of ln(h) of the macropore domain; θ r , h m , and σ are parameters of the matrix domain; and W is a constant in the range [0, 1).

Theoretical development of novel bimodal saturated hydraulic conductivity
We report on further adaptations to the physical model of Pollacco et al. (2013) to suit it to dual-porosity structured soils, which are common in New Zealand, solely from Kosugi hydraulic parameters describing θ (h).This involves rewriting the Romano bimodal θ (h) (Sect.3.1) and developing a novel bimodal K s model based on the modified bimodal θ (h) (Sect.3.2).

Modified Romano bimodal water retention curve
We propose a modified version of θ bim_rom (h) (Eq.9) that does not use the empirical parameter W .Our modified function, θ bim (h), is plotted in Fig. 1 and is computed as where θ s_mac is the saturated water content that theoretically differentiates macropore and matrix domains.
The shape of θ bim (h) is identical to that of θ bim_rom (h), but the advantage of θ bim (h) is that it uses the physical parameter θ s_mac instead of the empirical parameter W , and θ s_mac (≤ θ s ) is more easily parameterized than W particularly when there are no available data in the macropore domain.When we do not have data in the macropore domain, θ s_mac is determined by fitting the hydraulic parameters θ s_mac , θ r , h m , and σ of θ bim_mat (h) (Eq.11) solely in the matrix range (r < R mac or h > H mac ).Fig. 1 shows that R mac and θ s_mac delimit the matrix and the macropore domains and that r m of the Kosugi model is the inflection point of θ bim_mat (h) and r m_mac is the inflection point of θ bim_mac (h).

Novel bimodal saturated hydraulic conductivity model
Using θ bim (h), we propose a new bimodal K s_bim that is derived following K s_uni (Eq.7) but for which we add a macropore domain: ) dS e (11b) where r macropore is r ≥ R mac and r matrix is r < R mac .
The r matrix of Eq. ( 14) is derived from Eq. (2b): and r macropore is computed similarly as We introduced r matrix (Eq.16) and r macropore (Eq.17) into K s_bim (Eq.13), giving the equation for K s_bim : In Eq. ( 19), r m_mac is replaced by Y /h m_mac and r m is replaced by Y /h m .Note that the bimodal K s model requires that the flow in the macropore domain obeys the Buckingham-Darcy law.Therefore, this model's performance may be restricted in cases of non-Darcy flow, such as non-laminar and turbulent flow, which may occur in large macropores.
In this study, σ _mac is not derived from measured θ (h) because measured data in the macropore domain are not always available, and so it will be treated as a fitting parameter.As discussed above, θ s_mac , θ r , σ , and h m are optimized with θ uni (h) measurement points only in the matrix range (r < R mac or h > H mac ), which means that θ s is not included in the observation data.In summary, K s_bim requires optimization of the parameters τ 1 , τ 2 , τ 3 , and τ 1_mac , τ 2_mac , τ 3_mac , h m_mac , and σ _mac (if no data are available in the macropore domain).The theoretically feasible range of the parameters of K s_bim is shown in Table 3.
One of the limitations of the New Zealand data set is that it has no θ (h) data points in the macropore domain.The closest data point near saturation is θ (h = 50 cm), which is in the matrix pore space.Carrick et al. (2010) found that H mac ranges from 5 to 15 cm, with an average H mac = 10 cm, which corresponds to a circular pore radius of R mac = 0.0149 cm (e.g.Jarvis, 2007;Jarvis and Messing, 1995;Messing and Jarvis, 1993).Therefore, to reduce the number of optimized parameters we make the following assumption: where P m_mac is a fitting parameter greater than 1.We found the fitted value of P m_mac was 2.0; however, this fitted parameter was very broadly determined.The cause might be To avoid any unnecessary overlap of θ bim with θ bim_mat 1 > τ 1 > τ 1_mac ≥ 0 Flow in the macropore domain (larger pores) is expected to be straighter than in the matrix domain (smaller pores) due to reduced crookedness of the path 1 > τ 2 > τ 2_mac ≥ 0 It is expected that the shape of the "microscopic capillary tube" of the macropore domain (larger pores) is more perfectly cylindrical than in the matrix domain (smaller pores) The macropore domain has larger pores, and therefore it is assumed that the pores are better connected than the matrix pores that we are optimizing σ _mac , and therefore h m_mac and σ _mac might be linked.Linked parameters (Pollacco et al., 2008a(Pollacco et al., , b, 2009) ) mean that there is an infinite combination of sets of linked parameters h m_mac and σ _mac which produces values of objective function close to that obtained with the optimal parameter set and for which there exists a continuous relationship between h m_mac and σ _mac .Further research needs to determine if having more data in the macropore domain would reduce the cause of non-uniqueness.To illustrate h m_mac , the equivalent r m_mac point is shown in Fig. 1, where r m_mac is the inflection point of the macropore domain.Figure 1 also shows that the matrix and the macropore domains meet at R mac (H mac ).

Measurement of physical soil properties
The soil data used in this study were sourced from two data sets.In the first data set (Canterbury Regional Study; Table 4), soils were derived from eight soil series on the postglacial and glacial alluvial fan surfaces of the Canterbury Plains (Webb et al., 2000).The soils varied from shallow, well-drained silt loam soils to deep, poorly drained clay loam soils.The second data set was derived from the Soil Water Assessment and Measurement Programme to physically characterize key soils throughout New Zealand in the 1980s.Soils selected from this data set are listed by region in Table 4 and were selected from soils formed from sediments derived from indurated sandstone rocks, because this is the most common parent material for soils in New Zealand and has a reasonably representative number of soils analysed for physical properties.
The cores for particle size analysis and measurement of θ(h) had diameters which ranged from 5.5 to 10 cm diameter and height which varied from 5 to 6 cm.The 5 and 10 kPa measurements of the θ (h) were derived using the suction table method as per Dane and Topp (2002), following the NZ Soil Bureau laboratory method (Gradwell, 1972).The 20 to 1500 kPa of the θ (h) were measured using the pressure plate method as per Dane and Topp (2002), following the NZ Soil Bureau method (Gradwell, 1972).The laboratory analysis for particle size followed Gradwell (1972).
The total porosity, φ, described in Eq. ( 5) contains uncertainties from the measurement methods, where φ is derived from separate measurements of particle density and bulk density, rather than being directly measured.The uncertainty in φ measurements appeared to have reduced the demonstrated benefits of using K s_bim instead of K s_uni , which strongly relies on φα − θ s_mac and may have caused the optimal α to be 0.98 and not the commonly accepted value of 0.95 (Rogowski, 1971;Pollacco et al., 2013;Haverkamp et al., 2005;Leij et al., 2005).
The K s data used were collected and processed at a time when the best field practices in New Zealand were still being explored.K s was derived using constant-head Mariotte devices (1 cm head) from three to six cores (10 cm diameter and 7.5 cm thickness) for each horizon.The log 10 scale value of the standard error of the replicates of the measurements is shown in Fig. 2, which shows large uncertainty in the measurements (up to 3 orders of magnitude).This uncertainty is due to a. measurements of θ (h) and K s being taken on different cores, which caused some mismatch between θ (h) and K s , resulting in 16 outliers that negatively influenced the overall fit of the K s model having to be removed from the data set; b. side wall leakage of some cores, which led to K s values that were too high (Carrick, 2009), resulting in six samples with unusually high K s having to be removed from the data set; c. misreporting low K s since the measurements of K s were halted when conductivity was less than 0.1 cm day   2002;Anderson and Bouma, 1973) that were evident in dyed samples; we therefore removed measuredK s replicates that were too high and showed evidence of macropore abundance by having values of θ s − θ s_mac >0.05.
We therefore selected 235/262 samples (90 %) and removed only 27 outliers, which is minimal compared, for instance, to the UNSODA (Leij et al., 1999;Schaap and van Genuchten, 2006) and HYPRES databases (Wösten et al., 1998(Wösten et al., , 1999;;Lilly et al., 2008), which are used for the development of PTFs such as the ROSETTA PTF (Patil and Rajput, 2009;Rubio, 2008;Young, 2009), and which were found to contain a large number of outliers.Using these databases, Pollacco et al. (2013) selected only 73/318 soils (23 %), which complied with strict selection criteria prior to modelling.Note that the K s observations in the topsoils have greater variability than in the subsoil layers (Fig. 2).This is because topsoils are more disturbed by anthropogenic disturbance and biological activity.Therefore, the topsoils also have a greater abundance of macropores and therefore are more prone to error when the sampling is performed with a small core size that does not contain a representative volume of the macropore network.

Inverse modelling and goodness of fit
The parameterization of the model was performed in two consecutive steps: 1. Optimization of θ s_mac , θ r , h m , and σ of the unimodal Kosugi θ bim_mat (h) (Eq.11) was performed by matching observed and simulated θ (h) in the range h < H mac (as discussed, θ s is not included in the observation data since we did not have data in the macropore domain).
The feasible ranges of the Kosugi parameters are described in Table 1.
The inverse modelling was performed in MATLAB using AMALGAM, which is a robust global optimization algorithm (http://faculty.sites.uci.edu/jasper/sample/)(e.g.ter Braak and Vrugt, 2008).For each step, we minimized the objective functions described below.The objective function, OF θ , used to parameterize Kosugi's θ (h) at the following pressure points (5, 10, 20, 40, 50, 100, 1500 kPa), is described by where the subscripts sim and obs indicate simulated and observed values, respectively.P θ is the set of predicted parameters (θ s_mac , θ r , h m , σ ) and P ower is the power of the objective function.The computation of K s_bim requires θ (h) to be accurate near saturation, when the drainage is mostly from large pores, and to achieve this balance we found by trial and error that best results are achieved when P ower = 6.The parameters of K s_uni and K s_bim models were optimized by minimizing the following objective function OF ks : where the subscripts sim and obs indicate simulated and observed values, respectively.P ks is the vector of the unknown parameters.The log transformation of OF ks puts more emphasis on the lower K s and therefore reduces the bias towards larger conductivity (e.g.van Genuchten et al., 1991;Pollacco et al., 2011).Also, the log transformation considers that the uncertainty in measured unsaturated hydraulic conductivity increases as K(θ ) increases.
The goodness of fit between simulated (K s_uni or K s_bim ) and observed K s was computed by the RMSE log10 : where N is the number of data points.
The following transformation was necessary to scale the parameters to enable the global optimization to converge to a solution: where T 1 is a transformed tortuosity τ 1 .Introducing Eq. ( 19) into K s_bim Eq. ( 14) gives

Results and discussion
We report on (1) the suitability of the K s_uni model (developed with European and American data sets, Pollacco et al., 2013) to predict K s for New Zealand soils experiencing large uncertainties, as shown in Fig. 2; (2) improvements made by stratifying the data with texture and topsoil/subsoil; and (3) enhancements made by using the bimodal K s_bim instead of the unimodal K s_uni .

Improvement made by stratifying with texture and topsoil/subsoil
It was expected that stratifying with texture and topsoil/subsoil (layers) should improve the predictions of K s to only a modest degree.This is because K s_bim and K s_uni are physically based models that are based on pore-size distribution, and therefore stratifying with soil texture or topsoil/subsoil is not likely to provide extra information.For instance, Arya and Paris (1981) showed that there is a strong relationship between pore-size distribution and the particlesize distribution, and therefore adding soil texture information should not improve the model.As expected, no significant improvements were made by stratifying with soil texture compared with a model that groups all texture classes (loam and clay) and layers (topsoil and subsoil) (overall improvement of 3 %) (Table 5).However, a significant improvement was made by stratifying by layer (topsoil and subsoil) (overall improvement of 23 %), and therefore the remaining results are presented by stratifying by layer.These results are obtained because topsoils have higher macropores and a smaller tortuous path than that in subsoil, as demonstrated by τ 1_top > τ 1_sub or T 1_top < T 1_sub , τ 2_top > τ 2_sub , τ 3_top > τ 3_sub (Table 6).It is important to note that tortuosity decreases as τ gets closer to 1.

Improvement made by using K s_bim instead of K s_uni
Figure 3 shows an acceptable fit between K s_bim and K s_obs (RMSElog 10 = 0.450 cm day −1 ), recognizing that the observations contain large uncertainties since the measurements were taken by using small cores (Sect.4.1).The overall improvement made by using K s_bim is somewhat modest (5 % for all soils).As expected, the reasonable improvement is greater for topsoil containing higher macroporosity (12 % improvement) than for subsoil (4 % improvement) (Table 6).This is because topsoil has higher macropore θ mac (θ s − θ s_mac ) (Table 7) caused by earthworm channels, fissures, roots, and tillage than subsoil.The RMSElog 10 of K s_uni for subsoil is 0.47 cm day −1 (Table 6), which is slightly worse compared to the RMSElog 10 of 0.420 cm day −1 when using UNSODA and HYPRES data sets (Pollacco et al., 2013).
The reason K s_bim shows smaller-than-expected improvements compared to K s_uni requires further investigation and testing with a data set containing fewer uncertainties.One plausible explanation is that K s_bim is highly sensitive to θ s , computed from total porosity φ (Eq.6), which had inherent measurement uncertainties (Sect.4.1).In addition, the possible existence of non-Darcy flow in large biological pores may decrease the outperformance of the bimodal model over the unimodal model.

Optimal tortuosity parameters
The optimal tortuosity parameters of K s_bim and K s_uni (Table 6) show that the optimal parameters are within the physically feasible limits, except for τ 3_mac parameters of the subsoil, which are greater than τ 3 .This is understandable because Pollacco et al. (2013) found τ 3 not to be a very sensitive parameter.As expected, T 1_mac is smaller than T 1 (τ 1_mac > τ 1 ), which suggests that the tortuosity parameters have a physical meaning.
The estimated value of the unimodal T 1 parameter K s_uni derived from the UNSODA and HYPRES data sets (T 1 = 0.1) (Pollacco et al., 2013) is very different from the value estimated in this present study (T 1 = 6.5).Cichota et al. (2013) also reported that PTFs developed in Europe and the USA were not applicable to New Zealand.The reasons  7. Descriptive statistics of the optimized θ mac (θ s −θ s_mac ), θ s , h m , and σ Kosugi hydraulic parameters.The bar represents the average value, SD is the standard deviation, and N the number of measurement points.why these PTFs are not directly applicable to New Zealand require further investigation.

Uncertainty of the bimodal saturated hydraulic conductivity model predictions
The practical application of the bimodal saturated hydraulic conductivity model, K s_bim , to New Zealand soils requires a model for the uncertainty of the resultant predictions, since it is then possible to attach a value for the uncertainty of future predictions of K s .In a conventional parametric statistical model, the uncertainty model follows from the structure of the fitting model itself.In the present work, K s is estimated using an inverse model and this has no associated functional uncertainty model.For this reason, the uncertainty is derived empirically by fitting a relationship between the transformed residuals of the model (the log-transformed measured K s minus the log-transformed estimated K s ) as a function of the log-transformed estimated K s .Although the uncertainty model could be derived from all the soils in the study, this process results in a pooled estimate for uncertainty (e.g.aggregated root mean square error).However, it has been observed that topsoils and subsoils have different uncertainty behaviour for the estimated K s , so it is desirable to include an indicator variable to determine whether the soil is a topsoil or not.In explicit form, log 10 K s obs − log 10 K s sim = a 1 L + a 0 + , where a 0 and a 1 are fitting constants, L is an indicator variable specifying whether the soil is a topsoil (value of 1) or a subsoil (value of 0), and is the uncertainty distribution.
The distribution of the uncertainty could take a number of forms, but there is no obvious choice, except that one might expect the distribution central measure to be unbiased.To avoid an explicit distribution assumption, we fitted a conditional quantile model (Koenker, 2005) for the transformed residuals, based on the τ quantile, where τ = 0.5 corresponds to the conditional median, and τ = 0.025 and τ = 0.975 correspond, respectively, to the 2.5 and 97.5 % quantiles and thus together describe the 95 % containment interval of the residuals.
The conditional quantile model Eq. ( 25) was fitted using τ = 0.5, 0.025 and 0.975 (Table 8).The results suggest a strong dependence of the scale of the residuals on whether the soil is a topsoil or not, but the size of the 95 % residual containment interval is not dependent on the simulated K s .Notably, the confidence interval for the fitted median (τ = 0.5) quantile model suggests that the uncertainty distribution median is unbiased; thus, predictions from K s_bim show no propensity for bias, which is a desirable result.
Another way to illustrate the uncertainty model is to plot the observed log 10 K s_obs against the estimated logK s_bim , with the fitted median, lower, and upper 95 % quantile lines, as shown in Fig. 4. The width of the 95 % containment interval for the residuals is narrower (i.e. the predictions appear to be more accurate) for topsoils.The quantile estimates for the conditional median of both topsoil and subsoil are also shown in Fig. 4, with the shaded region showing the 95 % confidence interval of the median estimate.The shaded region covers the 1 : 1 line in Fig. 4, and thus there is no compelling evidence that the median residual distribution is biased.A key outcome of this research will be to provide direction for future field studies to quantify soil water movement attributes of New Zealand soils and to prioritize which measurements will have the greatest value to reduce the uncertainty in modelling of the soil moisture retention and hydraulic conductivity relationships.Recommendations are to evaluate the spatial representativeness of the current soil physics data set and undertake more measurements of hydraulic conductivity and soil water retention on key soils; use larger cores for measurements of hydraulic conductivity; take measurements of the moisture retention curve and saturated hydraulic conductivity on the same sample; provide more accurate measurements of total porosity; conduct near-saturation measurements of θ(h) and K(θ ) to better characterize the macropore domain, which is responsible for preferential flow behaviour; and make more accurate measurements on slowly permeable soils (< 1 cm day −1 ), which are important for management purposes but are not well represented in the current databases.

Conclusions
We report here on further adaptations to the saturated hydraulic conductivity unimodal model to suit it to dualporosity structured soils, by computing the soil water flux through a continuous function of a modified version of the Romano et al. (2011) θ (h) dual pore-size distribution.The shape of the Romano θ (h) distribution is identical to the modified θ (h), but the advantage of the developed bimodal θ (h) is that it is more easily parameterized when no data are available in the macropore domain.
The stratification of the data with texture only (loam or clay) slightly improved the predictions of the K s model, which is based on pore-size distribution.This gives us confidence that the K s model is accounting for the effect of these physical parameters on K s .A significant improvement was made by separating topsoils from subsoils.The improvements are higher for the topsoil, which has higher macroporosity caused by roots and tillage compared to subsoils.The reason why a model with no stratification is not sufficient is unclear and requires further investigation.
The improvements made by using the developed bimodal K s_bim (Eq.20) compared to the unimodal K s_uni (Eq.8) are modest overall, but, as expected, greater for topsoils having larger macroporosity.Nevertheless, an acceptable fit between K s_bim and K s_obs was obtained when due recognition was given to the high variability in the measured data.We expect K s_bim to provide greater improvement in K s predictions if more θ (h) measurements are made at tensions near saturation and if measurements are made on larger cores and with more accurate measurements of porosity.

Figure 2 .
Figure 2. Uncertainty of the standard error of the observed K s in topsoil and subsoil.The lines in the box show upper and lower quartiles, the median (red), and mean (green).Whiskers show values within 1.5 times the quartile spread; values outside this range are shown as plotted points.

KFigure 3 .
Figure3.Plot between K s_obs against K s_bim and K s_uni for topsoil and subsoil.The dotted line refers to the 1 : 1 line.

Figure 4 .
Figure 4. Error of K s_bim plotted against K s_obs for topsoil and subsoil.The solid line refers to the median line for each group, the dashed line refers to the upper or lower 95 % confidence interval lines, the dotted line refers to the 1 : 1 correspondence line, and the shaded region is the 95 % confidence interval of the median estimate.

Table 1 .
Feasible range of the Kosugi parameters and θ 5 (which is θ measured at 5 kPa).

Table 2 .
Description of the tortuosity parameters.

Table 3 .
Theoretical constraints of the K s_bim model.

Table 4 .
Soil series and classification.

Table 5 .
The RMSE log10 reported by using K s_bim and K s_uni models, by stratifying the data with/without texture and layers.

Table 8 .
Summary of the quantile regression fit of the logtransformed residuals.