Validation of Soil Survey Estimates of Saturated Hydraulic Conductivity in Major Soils of Puerto Rico

: Ranges or “classes” of probable saturated hydraulic conductivity values ( K sat ) are listed for all soil series in USDA-NRCS Soil Survey reports. Listed values are not measured, but rather estimated from other soil properties using a pedotransfer function (PTF). To validate the PTF, we compared estimated K sat classes with measured values in various horizons of nine major soil series of Puerto Rico. For each horizon, a minimum of 9 and usually 16 K sat measurements were made with Guelph permeameters near locations where soil pedons had been thoroughly described. In most horizons, K sat was log-normally distributed. The ratios of K sat values corresponding to one geometric standard deviation above and below the mean were usually less than 10, which is the ratio of upper and lower class boundaries in the K sat classiﬁcation system. For most horizons, measured K sat values were distributed among the rated K sat class and the next higher class, indicating that the PTF systematically underestimated the K sat distributions, but by less than an order of magnitude. From the point of view of soil and water management decisions requiring conservative K sat estimates, the PTF estimates appeared reasonably conservative without deviating from actual values so as to limit the usefulness of the estimates.


Introduction
One of the fundamental physical properties of soil is its saturated hydraulic conductivity (K sat ), defined as the rate of water movement through saturated soil under a unit hydraulic gradient. Knowledge of K sat is essential for predicting the magnitude of environmental processes such as water infiltration, runoff, and soil erosion, and for designing irrigation, drainage, and land-applied waste disposal systems [1][2][3][4][5]. Engineering properties of soil, such as consolidation rate, fluidization, piping, and embankment stability, are strongly influenced by the capacity of water to move through the soil, and, hence, by K sat [6,7]. For these reasons, information on K sat is a standard component of Soil Survey reports.
Field and laboratory determinations of K sat are expensive and time consuming, and require large numbers of measurements in order to account for high coefficients of variability which can exceed 400 percent [8]. The variability is compounded when K sat measurements are separated by large spatial distances [8][9][10], or taken at different times [11,12]. This has led to estimates of K sat based on pedotransfer functions (PTFs), which estimate values from correlated soil properties that are readily available in Soil Survey reports, such as texture, bulk density, organic matter, particle size distribution, and structure descriptions. Early evidence of the usefulness of PTFs for estimating K sat came from the work of O'Neal [13,14], who developed a set of field clues for estimating K sat . He found good agreement between measured and estimated K sat values in 68 percent of a total of 271 soil horizons examined.
To account for intrinsic soil variability, he assigned each of the soil horizons to one of seven "percolation classes" defined by ranges of probable K sat values, rather than assigning a single average value to each horizon. A similar "class" approach for estimating K sat in soils based on field textural and structural observations was adopted by McKeague et al. [15], who observed that 87 percent of estimated K sat values from 78 soil horizons deviated from measured values by one K sat class or less. Beginning with such seminal works, an enormous body of literature on PTF methodology has flourished [16][17][18][19][20][21]. These vary in complexity, from simple class systems to highly complex schemes based on neural networks.
One of the simplest yet most widely used PTFs for estimating K sat in the continental USA and its territories is described in Section 618 of the National Soil SurveyHandbook published by the USDA-NRCS [22]. The structure of the PTF is similar to that outlined by Pachepsky and Park [23]. Soils are separated into groups based on their location in the classical USDA textural triangle and bulk density groups defined broadly as "high", "medium" and "low". Depending on soil location in the textural triangle and bulk density grouping, soils are assigned to classes of probable K sat values, each varying over an order of magnitude. A list of auxiliary or overriding criteria, consisting primarily of soil structural descriptions, is also provided. Whenever one of the overriding criteria is encountered, a K sat class specific to that overriding criterion is assigned regardless of soil texture and bulk density. Estimates of K sat class based on this PTF are available for all soil series of the US, including Puerto Rico, and are listed in the NRCS National Soil Information System (NASIS) database. The different classes are defined in Table 1. Published studies are available comparing predictions of large numbers of PTF models to data sets of thousands of field and laboratory K sat measurements [17][18][19][20][21][22][23]. However, in surveying this literature, it is difficult to find studies evaluating the specific USDA-NRCS pedotransfer function described above. Pachepski and Park [23] compared a similar model based on texture and bulk density to several thousand K sat measurements in the USA. They found that the predictive accuracy of the PTF was not high, and yet was comparable with estimates obtained from far more detailed soil information using sophisticated machine learning methods.
Gupta et al. [21] assembled a global database of soil saturated hydraulic conductivity (SoilKsatDB) involving 13,267 K sat measurements from 1910 sites for use in PTF modeling. They commented that PTFs derived for soils in temperate soils may not be suitable for estimating K sat in tropical regions. This is of particular relevance to K sat estimates in Puerto Rico and other tropical islands of the Caribbean basin, which are based on PTF models developed primarily in temperate regions.
The objective of this research was to validate K sat estimates based on the USDA-NRCS pedotransfer function for major soils of Puerto Rico occupying large land areas and holding key positions in the soil classification system. The estimates are listed in NASIS databases for these soils. Validation was achieved by comparing estimated K sat classes with multiple measured values in the same soils.

Description of the Project Area
Saturated hydraulic conductivity was measured in situ in various horizons from 9 soil series of Puerto Rico, which had previously been thoroughly characterized by NRCS Soil Survey personnel. The soil series and their classifications according to the USDA Soil Taxonomy system are given in Table 2, together with the approximate Reference Soil Group in the World Reference Base (WRB) system [24]. The K sat measurement sites for each soil series were located as close as possible to pedons which had previously been characterized by Soil Survey personnel. The geographic locations and USDA-NRCS Pedon Identification Numbers of these pedons are given in Table 3, and a corresponding map is given in Figure 1. To minimize temporal effects on K sat , the measurements were made in pasture lands that had been undisturbed for several years prior to sampling. All measurements at a given site were made on the same day. Table 3. Coordinates of soil pedons nearest K sat measurement sites.

Experimental Techniques
At each site, the usual procedure was to set up a 4 × 4 sampling grid (for a total of n = 16 grid points), with a distance between grid points of approximately 6 m. The only exceptions were the sampling site for the Pandura series, where a 3 × 4 sampling grid was used, and the Bahia site, with a 3 × 3 sampling grid. A relatively close 6 m spacing between Ksat measurements at each site was used in order to capture intrinsic (random) soil variability near the pedon and minimize spatially correlated variations. At each grid point, Ksat was measured in situ using variants of the Guelph permeameter method [25,26], with the specific variant depending on whether Ksat was measured in surface or subsurface horizons.
For measurements in surface horizons, the top 2 cm of soil was removed to eliminate vegetation and debris. Sharpened PVC or stainless steel cylinders 10 cm long and of 10 cm interior diameter were driven into the ground to a depth (d) of 4 cm ( Figure 2). Soil immediately in contact with the inner and outer sides of the cylinder was pressed against the sides with a sharp stick to eliminate air gaps between the soil and cylinder walls. Ponded water at constant head (H) was then established inside the cylinders using a Guelph permeameter (Soilmoisture Equipment Corp., Santa Barbara, CA, USA), and the volumetric outflow rate (q) (volume time ⁄ ) from the permeameter was measured once steady state infiltration was reached, usually within 30-60 min.

Experimental Techniques
At each site, the usual procedure was to set up a 4 × 4 sampling grid (for a total of n = 16 grid points), with a distance between grid points of approximately 6 m. The only exceptions were the sampling site for the Pandura series, where a 3 × 4 sampling grid was used, and the Bahia site, with a 3 × 3 sampling grid. A relatively close 6 m spacing between K sat measurements at each site was used in order to capture intrinsic (random) soil variability near the pedon and minimize spatially correlated variations. At each grid point, K sat was measured in situ using variants of the Guelph permeameter method [25,26], with the specific variant depending on whether K sat was measured in surface or subsurface horizons.
For measurements in surface horizons, the top 2 cm of soil was removed to eliminate vegetation and debris. Sharpened PVC or stainless steel cylinders 10 cm long and of 10 cm interior diameter were driven into the ground to a depth (d) of 4 cm ( Figure 2). Soil immediately in contact with the inner and outer sides of the cylinder was pressed against the sides with a sharp stick to eliminate air gaps between the soil and cylinder walls. Ponded water at constant head (H) was then established inside the cylinders using a Guelph permeameter (Soilmoisture Equipment Corp., Santa Barbara, CA, USA), and the volumetric outflow rate (q) (volume/time) from the permeameter was measured once steady state infiltration was reached, usually within 30-60 min.

Experimental Techniques
At each site, the usual procedure was to set up a 4 × 4 sampling grid (for a total of n = 16 grid points), with a distance between grid points of approximately 6 m. The only exceptions were the sampling site for the Pandura series, where a 3 × 4 sampling grid was used, and the Bahia site, with a 3 × 3 sampling grid. A relatively close 6 m spacing between Ksat measurements at each site was used in order to capture intrinsic (random) soil variability near the pedon and minimize spatially correlated variations. At each grid point, Ksat was measured in situ using variants of the Guelph permeameter method [25,26], with the specific variant depending on whether Ksat was measured in surface or subsurface horizons.
For measurements in surface horizons, the top 2 cm of soil was removed to eliminate vegetation and debris. Sharpened PVC or stainless steel cylinders 10 cm long and of 10 cm interior diameter were driven into the ground to a depth (d) of 4 cm ( Figure 2). Soil immediately in contact with the inner and outer sides of the cylinder was pressed against the sides with a sharp stick to eliminate air gaps between the soil and cylinder walls. Ponded water at constant head (H) was then established inside the cylinders using a Guelph permeameter (Soilmoisture Equipment Corp., Santa Barbara, CA, USA), and the volumetric outflow rate (q) (volume time ⁄ ) from the permeameter was measured once steady state infiltration was reached, usually within 30-60 min.  Hydraulic conductivity K sat was calculated from the formula [16]: where the parameter G is a shape factor given by (2) and the cylinder configuration parameters a, d and H are defined in Figure 2. The parameter λ in Equation (1) with dimensions of length (cm), is the exponent in the hydraulic conductivity function where h (cm) is the suction head and K (h) is the soil hydraulic conductivity corresponding to h. The parameter λ has been used as a flow-weighted wetting front suction head in Green-Ampt infiltration models [25].
Estimates of λ for different texture/structural classes are described in Table 4. As all soils in this study were reasonably well structured, the parameter value λ = 8 cm was assumed. Although λ is estimated rather than directly measured, the resulting error in K sat estimates is usually small, within a factor of 2 of the actual value [25,26]. This uncertainty is small compared to within field variations of K sat , which can vary over 1 to 2 orders of magnitude. Table 4. Texture-Structure categories for visual estimation of λ. Adapted from Elrick et al. [25].

Texture-Structure Category λ * (cm)
Compacted, structureless, clayey, or silty materials such as landfill caps and liners, lacustrine or marine sediments, etc. 100 Most structured and medium textured materials; include structured clayey and loamy soils, as well as unstructured medium sands. This category is generally the most appropriate for agricultural soils.

25
Soils that are both fine textured and massive; include unstructured clayey and silty soils, as well as structureless sandy materials. 8 Coarse and gravelly sands; may also include some highly structured soils with large numerous cracks and biopores. 3 * The parameter λ in Table 9 and Equations (1), (3) and (4) is the inverse of another parameter (α) which was used by Reynolds et al. [25,26] in their description of Guelph permeameter theory. Our preference for using λ rather than α stems from the clear physical meaning of the former as a flow-weighted wetting front suction head in the Green-Ampt infiltration model [25].
For measuring K sat in subsurface horizons, a cylindrical hole of radius a = 5 cm was bored to the desired depth where K sat was to be measured. The bottom of the hole was flattened with a planing auger. To minimize soil smearing effects on K sat measurements, augering of very wet soils was avoided, and the walls of the auger hole were cleaned with a stiff brush which partially removed any smear layer. The permeameter tip was placed in the hole and water was allowed to infiltrate the soil until a steady state was reached under a constant water head (H) which was normally 5 cm, and 10 cm in very impermeable soils. The system is illustrated in Figure 3.
The hydraulic conductivity K sat was calculated from the measured steady state flux q as [26]: where C is a constant shape factor [25] that is dependent on λ and the ratio H/a. In our measurements, the well radius a was 3 cm and the hydraulic head H was either 5 or 10 cm.
where C is a constant shape factor [25] that is dependent on and the ratio / . In our measurements, the well radius a was 3 cm and the hydraulic head H was either 5 or 10 cm. The Guelph permeameter method was chosen as it is a theoretically well-founded field method that allows measuring Ksat in situ using relatively small amounts of water. A requirement is that a steady state water infiltration (q) must be achieved prior to measurement, but in most cases this is achieved within 30-60 min of infiltration. As illustrated in Figures 2 and 3, the infiltrating water is confined to a relatively small bulb-shaped wetting zone in the soil, which allows good depth resolution of Ksat measurements. The method can be applied using either one or two values of hydraulic head (H). In the single-head method, the infiltration rate (q) is measured at a given head (H) and the wetting front suction parameter is estimated from Table 4. This method typically gives a measurement accuracy of 0.5 ≤ ≤ 2 , where X is the measured value and Ksat is the true saturated conductivity value [24]. In the two-head method, the infiltration rate (q) is measured at two different infiltrometer heads (H), which allows the measurement of both Ksat and , and typically gives more accurate values of Ksat. However, this method consumes twice the amount of water and time as the single-head technique, and is mathematically ill-conditioned, which can cause unreasonable estimates of Ksat and [24,25]. Since our study involved large numbers of measurements to be completed in a single day at any given field site, requiring hand-carrying of all the necessary water, practical considerations dictated that the single-head method was preferable.

Results
Cumulative distributions of measured Ksat values in the different soil horizons are shown in Figure 4. Shapiro-Wilks tests on the distributions showed that they were nearly all log-normally distributed. Consequently, all subsequent parametric statistical tests were performed on log-transformed Ksat values. The Guelph permeameter method was chosen as it is a theoretically well-founded field method that allows measuring K sat in situ using relatively small amounts of water. A requirement is that a steady state water infiltration (q) must be achieved prior to measurement, but in most cases this is achieved within 30-60 min of infiltration. As illustrated in Figures 2 and 3, the infiltrating water is confined to a relatively small bulb-shaped wetting zone in the soil, which allows good depth resolution of K sat measurements. The method can be applied using either one or two values of hydraulic head (H). In the single-head method, the infiltration rate (q) is measured at a given head (H) and the wetting front suction parameter λ is estimated from Table 4. This method typically gives a measurement accuracy of 0.5K sat ≤ X ≤ 2K sat , where X is the measured value and K sat is the true saturated conductivity value [24]. In the two-head method, the infiltration rate (q) is measured at two different infiltrometer heads (H), which allows the measurement of both K sat and λ, and typically gives more accurate values of K sat . However, this method consumes twice the amount of water and time as the single-head technique, and is mathematically ill-conditioned, which can cause unreasonable estimates of K sat and λ [24,25]. Since our study involved large numbers of measurements to be completed in a single day at any given field site, requiring hand-carrying of all the necessary water, practical considerations dictated that the single-head method was preferable.

Results
Cumulative distributions of measured K sat values in the different soil horizons are shown in Figure 4. Shapiro-Wilks tests on the distributions showed that they were nearly all log-normally distributed. Consequently, all subsequent parametric statistical tests were performed on log-transformed K sat values.
For such distributions, the arithmetic mean log K from a sample of n = 1, 2, . . . , N values of log K sat values is given by where K geo is the geometric mean defined by The standard deviation of log K sat values, σ log K is defined by where σ geo is the geometric standard deviation, defined by where is the geometric mean defined by The standard deviation of log Ksat values, is defined by where is the geometric standard deviation, defined by   The parameters K 16.5 , K 50 and K 83.5 are K sat values corresponding to the cumulative distribution percentiles 16.5, 50, and 83.5, respectively. The parameters K 83.5 and K 16.5 constitute the upper and lower bounds of the middle 67 percent of all K sat values, i.e., the percent of values residing within one geometric standard deviation of the geometric mean. Table 5 lists the soil horizons evaluated, the number N of K sat measurements per horizon, the rated K sat class in the USDA-NRCS system, the measured geometric mean of K sat , the geometric standard deviation, the respective upper and lower geometric standard deviation bounds K 83.5 and K 16.5 , and the ratios K 83.5 K 16.5 . Statistical differences among the means of the distributions are indicated by the lower case letters in the second column. Soil horizon names followed by a given letter in common are not statistically different at the p < 0.05 level. Statistical differences were determined by performing ANOVA and Tukey tests on the log transformed K sat values. Each sampled soil horizon was considered a different experimental treatment, and the number of K sat measurements for that horizon were taken as the corresponding number of replications. Results show that 9 of the 14 soil horizons are not statistically different from one another. The two sandy Bahía soil horizons had K sat values significantly higher than all other soil series except Coto, a clayey Oxisol. The two soils with statistically lower K sat values than all others were Descalabrado and Fraternidad, both characterized by abundance of 2:1 clay minerals.

Fraction of sample measurements with lower Ksat values
Inspection of the right hand columns of Table 5 shows that the ratios of K sat values corresponding to one geometric standard deviation above and below the mean, i.e., K 83.5 /K 16.5 , were usually less than 10, which is the ratio of upper and lower class boundaries in the K sat classification system. This means that, in most cases, the assigned class bandwidth of an order of magnitude was sufficient to capture most of the variability of K sat values in a given sample. The two most important exceptions occurred in the case of the two soils with lowest K sat values, Fraternidad and Descalabrado, characterized by abundance of 2:1 clay minerals. In this case, the ratio K 83.5 /K 16.5 was approximately 20, indicating a considerably broader distribution than that assumed in the rating system. Inspection of Table 5 and Figure 5 shows a reasonably strong inverse relation between the dispersion parameter K 83.5 /K 16.5 and the geometric mean K geo .
However, although the measured K sat distributions were of comparable dispersion as the rated class boundaries, the measured K sat distributions were always shifted upwards relative to the rated K sat class, indicating that the PTF systematically underpredicted the actual values. This can be seen in Table 6, which gives the fraction of measured K sat values occurring within the expected (rated) class, and the fractions occurring in classes one or two orders of magnitude higher or lower than the expected class. For all soils, the bulk of measured K sat values occurred either in the rated class or in classes above it, with almost no measured values occurring below the rated class. For most soils, 90 percent or more of the K sat data were distributed among the rated class and the class immediately above it, indicating that the PTF underpredicted the actual K sat values by no more than one class or one order of magnitude. A significant exception occurred in the case of the two Humatas horizons, for which the PTF placed the rated class (0.1-1 µm s −1 ) one or two orders of magnitude below the actual K sat values. However, although the measured Ksat distributions were of comparable dispersion as the rated class boundaries, the measured Ksat distributions were always shifted upwards relative to the rated Ksat class, indicating that the PTF systematically underpredicted the actual values. This can be seen in Table 6, which gives the fraction of measured Ksat values occurring within the expected (rated) class, and the fractions occurring in classes one or two orders of magnitude higher or lower than the expected class. For all soils, the bulk of measured Ksat values occurred either in the rated class or in classes above it, with almost no measured values occurring below the rated class. For most soils, 90 percent or more of the Ksat data were distributed among the rated class and the class immediately above it, indicating that the PTF underpredicted the actual Ksat values by no more than one class or one order of magnitude. A significant exception occurred in the case of the two Humatas horizons, for which the PTF placed the rated class (0.1-1 μm s ) one or two orders of magnitude below the actual Ksat values.

Discussion
The above results confirm the need to assign a wide range in K sat values to a given class, as is currently recognized in the USDA-NRCS pedotransfer function. The measured dispersion of K sat values around the geometric mean was usually comparable to the 10fold dispersion implicit in the pedotransfer function, except in the case of soils with very low K sat values (soils with abundance of expandible 2:1 clay minerals) where a broader dispersion of values was observed. Other authors [27,28] have encountered spreads in K sat values commonly within ranges of one or two orders of magnitude.
For the most part, the USDA-NRCS pedotransfer function tended to underestimate K sat values, as indicated by the fact that most experimental values occurred either in the rated K sat class or in the class immediately above it. This result is similar to that of Sobieraj et al. [17], who compared K sat estimates of nine PTF models to measured values in a tropical watershed, and noticed that in the 0-0.1 soil depth interval, all PTF models slightly underestimated the K sat values. In our study, the difference between estimated and measured K sat values was usually no greater than one K sat class, comparable to the results of McKeague et al. [15] in soils from Canada and northeastern USA, and findings in the pioneering study by O'Neil [13,14] in soils from the USA.
This study only considered spatial variation of K sat measurements over short distances. To minimize the possibility of significant temporal effects, all measurements were performed in soil horizons that were under pasturewhich, other than grazing, had not been disturbed for at least three or four years. Kargas et al. [28] noted from their own work and other cited research that temporal variability of K sat can be small in soils covered by the same vegetation throughout the year and subject to no recent human intervention.
From the point of view of soil and water management decisions based on estimated K sat classes, it is often safer to underestimate K sat than to overestimate it. For example, in the case of land application of waste water or liquid manure, an underestimate of K sat will ensure that the recommended rate of source application (based on K sat ) will not overload the infiltration capacity of the soil, thereby reducing the risk of runoff and negative environmental impact. On the other hand, it is desirable that the estimated K sat class not be too far below the actual class, which would lead to an economically underdesigned system. For most soil horizon in this study, the USDA-NRCS estimation system appeared to satisfy these criteria by systematically underpredicting actual K sat values by a margin not exceeding one order of magnitude. The single exception to this observation was the Humatas soil, for which the rating system underpredicted K sat by one or two orders of magnitude.
We believe that the rather severe underprediction of K sat for the Humatas series, a dominant Ultisol in Puerto Rico, was probably due to over-weighting of texture and bulk density effects in the PTF, and not enough weighting of micro-structure effects associated with highly weathered clay minerals in Ultisols. Our suggestion is that for all Oxisols and Ultisols, which are highly weathered by definition, a K sat class of at least 1-10 µm s −1 should be assigned by default. This would be consistent with our results presented here, and also with a classic study by Lugo-Lopez et al. [29], where eighth hour infiltration rates measured with ring infiltrometers in Oxisols and Ultisols of Puerto Rico always exceeded 0.8 in hr −1 or 5.7 µm s −1 .