Modeling soil hydraulic properties for a wide range of soil conditions

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

available at www.sciencedirect.com

journal homepage: www.elsevier.com/locate/ecolmodel

Modeling soil hydraulic properties for a wide range of soil conditions Vincent Balland, Joseph A.P. Pollacco, Paul A. Arp ∗ Faculty of Forestry and Environmental Management, University of New Brunswick, 28 Dineen Drive, Fredericton, NB, Canada E3B 6C2

a r t i c l e

i n f o

a b s t r a c t

Article history:

This paper presents equations for estimating soil hydraulic properties such as bulk density

Received 15 March 2007

(Db), field capacity (FC), permanent wilting point (PWP), and saturated hydraulic conduc-

Received in revised form

tivity (Ksat ) from soil composition (sand, silt, clay, organic matter, coarse fragments), and

4 April 2008

with increasing soil depth, for a wide range of natural soil types and conditions. The equa-

Accepted 4 April 2008

tions, derived from New Brunswick and Nova Scotia soil survey reports, were evaluated and

Published on line 26 August 2008

extended to international databases, e.g., ISRIC (Db, FC, PWP) and UNSODA (Ksat ). Equation derivation was guided by theoretical considerations and regression techniques, to generate

Keywords:

realistic Db, FC, PWP and Ksat values with sand, clay, organic matter content and soil depth as

Field capacity

the only predictor variables, each with low inter-coefficient correlations. The results ensure

Permanent wilting point

that Db ≤ Dp (soil particle density), FC ≤ SP(soil saturation point), PWP ≤ FC, and Db, FC, PWP,

Saturated hydraulic conductivity

Ksat > 0. The best-fitted Db, FC, PWP, log10 (Ksat ), data variations are represented as follows:

Soil bulk density

R2 = 0.83, 0.96, 0.65 and 0.79, respectively. In addition, the associated root-mean-square-

Available water

residuals are similar to published values achieved with other approaches. The equations

Organic matter

can be used to quantify expected changes in available soil water, Db, FC, PWP, log10 (Ksat ),

Depth

and with varying organic matter content and texture. Examining the ISRIC data revealed

Texture

systematic coefficient variations for the Db equation by soil order.

Coarse fragments

1.

Introduction

Many soil management and engineering applications in forestry, agriculture, terrestrial ecosystem management and land reclamation require hydrologically feasible estimates for a number of soil properties that affect the retention and movement of water and dissolved substances through soils, for a wide range of conditions, from one extreme to another, i.e., from compacted to non-compacted soils, from organic to mineral soils, and for all organic matter and textural combinations, including coarse fragments, and with increasing soil depth. The properties that affect the hydrological behaviour of soils refer to soil pore space or the soil saturation point (SP), soil bulk density (Db), field capacity (FC), permanent wilting

∗

Corresponding author. E-mail address: [email protected] (P.A. Arp). 0304-3800/$ – see front matter © 2008 Elsevier B.V. All rights reserved. doi:10.1016/j.ecolmodel.2008.07.009

© 2008 Elsevier B.V. All rights reserved.

point (PWP), and saturated hydraulic conductivity (Ksat ,). Specific land-use applications generally build on estimating these properties for the mapping, projecting and executing best management operations on a daily year-round basis, whether this refers to, e.g., irrigation and water retention in crop management (agriculture, horticulture, forestry), soil water and groundwater engineering and management, soil conservation and rehabilitation, assessment of depth of heat and frost penetration into soils, and soil trafficability. From climate change and land-use change perspectives, hydrological soil properties are, for example, relevant to the forecasting of soil C turn-over rates and retention (Futter et al., 2007; Zhang et al., this issue). In turn, any changes in soil C may affect the soil hydrological properties due to the influence of plant litter and soil C

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

on microbial activities, and, hence, on soil structure and soil density. Even without these complications, mapping and modelling of Db, FC, PWP and Ksat reliably is problematic because these properties do not vary independently of one another, and vary strongly with changing soil texture (i.e., sand, silt and clay content), coarse fragment content (CF), soil depth, organic matter content (OM) and soil structure. For example, increasing soil OM levels usually increase the soil pore space and soil aggregation, and therefore affect Db, SP, FC, PWP, and Ksat (Hillel, 1982; Nemes et al., 2003; Rawls et al., 2003). Opposite changes can be expected with traffic-induced soil compaction. Any equations that address such changes need to ensure that the solutions do not move out of the expected feasibility ranges, so that FC < SP, PWP < FC, and Db < Dp, always. Equations that would produce non-feasible results jeopardize the numerical stability of hydrological models, especially when applied across widely varying soil conditions outside the immediate calibration range of these equations. In the literature, many equations for estimating SP, FC, PWP, and Ksat from soil survey data such as soil texture and organic matter content already exist. These equations are known as “point pedotransfer functions”, or point PTFs for short (Bouma, 1989; Schaap, 2006). In contrast, equations to estimate the parameters of water retention and conductivity in relation to soil matric potential (or head pressures) are known as parametric PTFs (Schaap, 2006). In general, PTF derivations involve, e.g., simple to multiple regression techniques, regression tree modeling, and artificial neural network (ANN) algorithms, coupled with varying optimization routines and other specialized software (Pachepsky et al., 1996; Schaap and Bouten, 1996; Tamari et al., 1996; Schaap and Leij, 1998; Pachepsky and Rawls, 2004; Schaap, 2006). In this, much attention is given to PTF accuracy, reliability, and cross-validations (Pachepsky and Rawls, 1999; Wösten et al., 2001; Gijsman et al., 2003). The continued compilation of regional and global soil databases is stimulating further PTF developments and investigations, from saturated to unsaturated soil conditions (Børgesen et al., 2007). In most cases, existing PTF derivations are empirical. The resulting functions are therefore limited to the range of the data from which the PTFs were generated. Outside these constraints, PTF projections may not only become inaccurate, but may also become infeasible. For example, the PTFs in Ritchie et al. (1987) produced unrealistic results for SAND = 0% or CLAY = 0%. The method of Saxton et al. (1986) was constrained to apply for soils with 5% ≤ SAND ≤ 30% and 5% ≤ CLAY ≤ 60%, in keeping with the available data range. Similarly, the PTFs generated by Rawls and Brakensiek (1989) and Rawls et al. (1991) were limited to soils with 5% ≤ CLAY ≤ 60% and 5% ≤ SAND ≤ 70%. Within these restrictions, there has also been a tendency to proliferate the number of predictor variables and coefficients. However, Schaap et al. (2004) have already shown that parsimonious PTFs tend to perform just as well as elaborate PTFs with many variables and coefficients. The primary objective of this paper is to present a systematic and generalizable framework for deriving physically constrained equations for Db, FC, PWP, and Ksat , from known texture and OM values in relation to soil depth for a wide range of natural soil conditions, as documented in soil survey reports and soil databases. The approach taken ensures

301

that the empirically derived equations would generate physically feasible solutions from fairly loose to fully compacted soils, from mineral soil with no organic matter to organic soils, from sandy soils to clay soils, etc. The equation derivation follows linear and nonlinear regression techniques with the added aim to retain only the most significant variables that are also accompanied with as precisely defined regression coefficients as possible. Secondary objectives deal with (i) comparing the precision of the approach with other recent PTF developments, (ii) applying the approach to the examination of several soil databases across various soil types, and (iii) presenting an application of the results to the forecasting of changes in soil density, depth, water retention, and permeability in response to soil organic matter amendment and subsequent soil re-compression, as affected by 0–100% clay content. The derivation of the proposed equations was at first restricted to soil data selected from two eastern Canadian Provinces, namely New Brunswick (NB) and Nova Scotia (NS), as found within the database of the Canadian Soil Information System (CANSIS). The general applicability of the equations was then evaluated with two global soil databases, namely the Universal Soil Database (UNSODA, Leij et al., 1996), and the database of the International Soil Reference and Information Centre (ISRIC, Tempel et al., 1996). In all of this, the soil. or – more specifically – any particular soil layer within a soil, is considered to represent an isotropic multi-phase system containing air, water, OM, mineral components finer than 2 mm (sand, silt and clay), and coarse fragments larger than 2 mm. Within-layer complications due to anisotropies resulting into flow fingering and channelling are therefore not addressed.

2.

Methods

2.1.

Data and metadata

Fig. 1 shows a schematic of the symbols used to describe and quantify the volume/weight relationships between Db, FC, PWP, and Ksat and the sand, silt, clay, OM and CF continuum from one extreme to the other. In this context, the “soil” subscript refers to the space between coarse fragments and roots, provided that space contains humified organic matter or minerals particles less than 2 mm. This definition is generally consistent with standardised soil sampling techniques, where soil cores are taken from the field, where the fine-earth fraction of these cores is separated from coarse fragments and roots, and where the sieved fraction is then analysed for texture, OM and other soil chemical and physical properties. Such cores may also be used for determining Db or SP, FC, PWP, and Ksat , with or without CF corrections (Soil Survey Laboratory Methods Manual, 2004). In the following, the terms “minerals” and “particles” are used to refer to the sum of SAND + SILT + CLAY = 1, and to the sum of OM + SAND + SILT + CLAY = 1, respectively. The weight fractions of materials that form the solid part of the soil are denoted by subscript “W”. For further details and definitions, and dealing with coarse fragments, see Appendix A. Data and metadata about Db, SP, FC, PWP, Ksat , soil texture, OM content, CF content, and soil depth were obtained

302

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

Fig. 1 – Overview of soil volume/weight symbols as used in the text: “V” and “W” represent volume and weight, respectively; “soil” refers to the fine earth fraction (any combination of sand, silt and clay minerals + soil organic matter + pores), but excludes coarse fragments (CF); “soil particles” refer to fine earth minerals and organic matter. Solid refers to fine earth + CF. Total volume refers to fine earth, CF and pores.

from NB and NS soil survey reports (CANSIS, 2000; Table 1), and from UNSODA and ISRIC data (Table 2). The NB and NS soil data were generated from undisturbed soils, and measurements were taken in situ, mostly referring to pristine forest conditions in a temperate to near-boreal climate, and mostly in a glaciated landscape with highly variable texture, OM, and drainage conditions. The ISRIC database represents a wide range of soil conditions in the USA. We considered that the spodosol and histisol entries in this database would in general be similar to the NB and NS soil data. We also used the entire ISRIC database for analytical purposes. The global UNSODA database was used to provide data for checking the general applicability of the PTF formulation for Ksat , as generated from the NB and NS data (the ISRIC database lacked this information). In some of the NB and NS references (Table 1), it was stated whether the soil samples had been disturbed as part of the sampling protocol for determining Db, FC, PWP, and Ksat (1), or were left undisturbed for the same purpose (2).

Table 1 – New Brunswick and Nova Scotia soils data chosen for developing generalizable Db, FC, PWP and Kset expressions Data sources Bulk density Field capacity Permanent wilting point Saturated hydraulic conductivity

a

(1) and (2) (1)a and (2) (1)a and (2) (2)–(5)

Number of data 237 227 209 112

(1) Langmaid et al. (1976); (2) Rees et al. (1992); (3) Mellerowicz et al. (1993); (4) Fahmy and Rees (1989); (5) Webb et al. (1991). a

For (1), it was decided to use data corresponding to undisturbed soil samples only. We compared bulk density measured in the field to bulk density calculated from pore space as determined in the laboratory and from the particle density, via Eq. (2). Only data with reasonable agreement (±0.1 g cm−3 ) between recorded and calculated field bulk densities were accepted.

In other references, field sampling and soil analytical procedures were not mentioned ((3)–(5)). For this reason, some of the data were either used or not, depending on whether decisions regarding the true nature of the data could be ascertained. Knowing this is particularly relevant for evaluating existing soil bulk density values, because these densities can be expected to depend on soil conditions at the time of sampling: while soil mass remains fixed, soil volume and soil bulk density may change due to compression pressures, ice conditions, and water content (Blake and Hartge, 1986; Pachepsky et al., 2001; Heuscher et al., 2005). For the purpose of this study, we use Db as reported for the field condition, or at or near field capacity, and apply the resulting values for estimating SP, FC, PWP and Ksat in an unmodified manner.

2.2.

Units

The units of the water content data varied among these references. For example (Table 1), water contents at saturation, field capacity and permanent wilting point were reported as “volume%” ((3) and (5)) or “oven-dried weight%” ((1) and (2)). In one report, this important qualifier was not mentioned (4). For the analysis of this study, all data for soil moisture, sand, silt, clay, and organic matter (NB, NS, ISRIC, UNSODA) were expressed in terms of “oven-dried weight, g g−1 , weight%”. All density units refer to g cm−3 , and Ksat units refer to cm h−1 unless otherwise stated.

2.3.

Analytical procedures

Regression analyses were done using the linear and nonlinear options of SPSS to determine how Db, FC, PWP, log10 (Ksat ) could be best predicted from SAND, SILT and CLAY content, OM, CF, and soil depth. Only those variables which would conform with general expectations and would also give the best overall fit with the proposed model formulizations were retained, e.g., high sand content should have a strong and

Table 2 – Summary for soil data used for developing generalized Db, FCWsoil , PWPWsoil and Ksat expressions Data source

Depth (cm)

Db (g cm−3 ) Sand

31.8 32.4 370 0 129.5

ISRIC spodosols, histosols (n = 234) Mean 57.42 S.D. 43.35 Minimum 2.50 Maximum 211.00 ISRIC all data Mean S.D. n Minimum Maximum UNSODA (n = 459) Mean S.D. Minimum Maximum

68.60 58.60 13103 0.00 790

– – – –

Silt

Clay

OM 0.13 0.3 409 0.001 1

0.21 0.21 277 0 0.78

1.26 0.51 407 0.1 2.26

0.39 0.23 409 0 0.97

0.31 0.17 409 0 0.70

0.19 0.13 409 0 0.72

1.43 0.35 0.64 2.24

0.60 0.18 0.05 0.98

0.32 0.16 0.00 0.90

0.05 0.03 0.00 0.24

0.03 0.03 0.00 0.23

1.57 0.25 13103 0.30 2.24

0.36 0.25 13103 0.00 0.99

0.37 0.19 13103 0.00 0.91

0.25 0.17 13103 0.00 0.94

1.49 0.21 0.83 1.98

0.55 0.29 0.00 1.00

0.30 0.22 0.00 0.81

0.15 0.13 0.00 0.63

CF

SP

FC

PWP

0.84 1.28 370 0.12 6.9

0.33 0.31 396 0.02 2.1

0.104 0.104 306 0.01 1.12

– – – –

0.40 0.25 0.07 1.46

0.22 0.16 0.03 0.99

0.065 0.05 0.01 0.29

0.01 0.02 13103 0.00 0.53

– – – –

0.35 0.17 13103 0.03 5.29

0.24 0.13 13103 0.01 3.04

1.42 1.51 0.00 7.31

– – – –

0.43 0.08 0.21 0.68

– – – –

0.13 0.07 12872 0.00 0.68

– – – –

13.9 25.6 123 0.01 157

– – – –

– – – –

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

NB, NS data Mean S.D. n Minimum Maximum

Ksat (cm h−1 )

Weight fractions

15.0 43.1 0.00 437.8

303

304

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

positive influence on saturated hydraulic conductivity. The resulting linear equations were then modified to ensure that the final equations comply with expected Db, FC, PWP and Ksat values and ranges for soils with very low density such as peat and organic soils, for mineral soils with very high density such as compacted clays, and for sandy soils with their compaction restrictions. The equations were required to produce feasible solutions only, i.e., all estimates for Db, FC, PWP and Ksat would have to be >0. Also, it was required that Db ≤ Dp, FC ≤ SP, and PWP ≤ FC without introducing discontinuities from one extreme to the other, i.e., any small change in texture, OM, CF, and soil depth should always produce a fully differentiable change in Db, FC, PWP, and Ksat . The statistical requirements that guided the equation derivations were as follows:

Dbsoil = 1.06 + 0.009 × DEPTH − 1.30 × OMWsoil

(1) However, when DEPTH is near zero, and OM = 1, this formulation generates Dbsoil = −0.24 g cm−3 . To be realistic: Db should always be greater than 0; for mineral soils, Dbsoil should be smaller than Dpsoil ; and for organic soils, Dbsoil should also be much smaller than 1.3 g cm−3 . The following nonlinear model keeps Dbsoil between 0 and Dpsoil , and (in comparison with the linear model) also provides a significantly better fit (R2 = 0.825): 1.23 + (Dpsoil − 1.23 − 0.75 × SANDWsoil ) Dbsoil =

• The number of predictor variables and adjustable coefficients would be reduced to the lowest possible number that would still give the highest correspondence (R2 values) between the actual and best-fitted values. • The equations would – as much as possible – produce evenly and normally distributed scatter plots of actual versus best-fitted values, with low root mean square errors (RMSE) and low bias (i.e., near-zero mean absolute errors (MAE). • The equations that would show high correlations among the best-fitted regression coefficients would be rejected, to reduce the numerical uncertainties associated with these coefficients as much as possible. • The derivation of equations followed these steps: ◦ first determine the overall soil particle density (Dpsoil ) based on mineral and OM contents; ◦ next, compute the soil bulk density, Db; ◦ then, compute the soil pore space SP next; ◦ thereafter, determine FC in relation to SP such that FC ≤ SP; ◦ also compute PWP in relation to FC such that PWP ≤ FC; ◦ finally, calculate Ksat or log10 (Ksat ).

In doing this, we used the NB and NS soil entries as a guide for deriving the linear and non-linear expressions for Db, FC, PWP and log10 (Ksat ). We then used the ISRIC subset for spodosols and histisols to check the resulting FC and PWP formulations, and the UNSODA data to check the resulting expression for log10 (Ksat ). Finally, we used the whole ISRIC database to generalize the equations for Db, FC, and PWP as much as possible.

Results and discussion

3.1.

Soil bulk density

Under natural conditions, field-state soil Db can generally be expected to increase with increasing soil depth, and to decrease with increasing OM. The best-fitted linear predictor model for the NS and NB soil Db data complied with these expectations, i.e.,

× (1 − exp(−0.0106 × DEPTH)) 1 + 6.83 × OMWsoil

(2)

Shown in Fig. 2 is the associated scatter plot of the bulk density data versus their best-fitted values. Also shown in Fig. 2 (left) are three-dimensional (3D) visualisations of how soil bulk density is therefore calculated to vary with changing sand and OMWsoil content, and soil depth. Several features are noteworthy: • The NB and NS data conform to the best-fitted calculations reasonably well with an RMSE of ±0.14 bulk density units. For the estimates concerning pore space or water content at field saturation by volume, this translates into an RSME precision of ±0.05 cm3 cm−3 (=0.14/2.65), or ±5%. This number is within the 3.3–6.3% RMSE range reported by Børgesen and Schaap (2005) for the best-fitted pressure-plate soil moisture content residuals at −1 kPa, generated with eight point PTFs or nine parametric PTFs (n = 1608, data set #2). The lower portion of this range was produced with 11–14 predictor variables involving 7 particle size determinations, 1–4 specifications for water content at given pressure heads. The upper portion was produced with 3 to 6 predictor variables involving sand, silt and clay, Db, OM, and a topsoil versus subsoil specification. Similarly, the analysis performed by Minasny et al. (2004) produced an RMSE range from 3.4 to 4.7% for water retention at any of seven measured water heads from 0 to 600 cm (number of soil samples = 310), and using 3–6 predictor variables. • Bulk density of soils with little or no organic matter should generally increases with depth, but less so for sandy soils than for silt and clay soils. Note that, with increasing soil depth, the formulation leads to: Dbsoil ≈

3.

(R2 = 0.733)

Dpsoil − 0.75 × SANDWsoil 1 + 6.83 × OMWsoil

(3)

This, in turn, implies that Dbsoil = Dpsoil when SANDWsoil and OMWsoil = 0, or when the combined silt and clay fraction = 1. This further implies that silt and clay soils would become fully compacted with increasing soil depth. • At the surface of the soil, the above equation reduces to Dbsoil ≈

1.23 1 + 6.83 × OMWsoil

(4)

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

305

Fig. 2 – Left: bulk density (Dbsoil ) in g cm−3 , according to Eq. (2). Top: predicted Dbsoil versus actual data. Center and bottom: 3D visualisations of how Dbsoil varies with increasing sand fraction, OM fraction, and soil depth (in cm). In this, sand fraction + silt fraction + clay fraction + OM fraction = 1. Right: field capacity (FCWsoil ) by weight, in %, according to Eq. (7). Top: predicted FCWsoil versus actual data. Center and bottom: 3D visualisations of how FCWsoil varies with increasing pore space (SPWsoil ), OMWsoil fraction and sand fraction.

This suggests that bulk densities at the surface of fine to coarse textured soils have, on average, a density of 1.23 g cm−3 when OMWsoil = 0, and that this density would decrease with increasing OMWsoil . • For organic soils with little mineral content, Eq. (2) becomes Dbsoil ≈

1.23 + (Dpsoil − 1.23)(1 − exp(−0.0106 × DEPTH)) 1 + 6.83 × OMWsoil

(5)

and predicts a Dbsoil of about 0.14 g cm−3 for purely organic soils such as peat, i.e., similar to values reported by Yu et al.

(2003), who reported that bulk density in a peat bog beneath the water table remained fairly constant between 0.1 and 0.15 g cm−3 (oven-dried basis, ±RMSE), to a depth of about 4, but rising sharply at that point to 0.2 g cm−3 , presumably on account of coprogenic inclusions. Above the water table, peat tends to self-compress (Schlotzhauer and Price, 1999), as also noted for compost piles (Agnew et al., 2003).

In general, bulk density profiles in soils are not only affected by soil depth, texture, and organic matter, but also

306

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

by recent and past soil compaction processes, such as glacial pressures, cycles of wetting and drying, freezing and thawing, surface traffic, and land-use practices. In the light of this, the above expressions for Db are intended to apply to soil bulk density profiles under natural conditions only, and for soils similar to those of the originating data in particular (i.e., Table 1). This equation, therefore, cannot be expected to apply to any engineered surfaces, where Db can and should be higher (in case of roads, trails, tracks, etc.) or lower (in newly established seed beds in gardens and on farm fields).

3.2.

Field capacity

The best-fitted linear predictor equation for the FCWsoil values of the NB and NS soils was FCWsoil = −0.33 + 0.32 × SANDWsoil − 1.53 × OMWsoil +1.0 × SPWsoil

(R2 = 0.983)

(6)

±0.104 cm3 cm−3 ; adding Db to these specifications reduced the RMSE to ±0.087 cm3 cm−3 . In another study, Pachepsky and Rawls (1999) reported an average RMSE of ±0.035 cm3 cm−3 and an RMSE range from about 0.018 to 0.066 cm3 cm−3 for soil moisture retention at −33 kPa for a wide range of soil groups (n = 447, but no organic soils), using linear and quadratic terms of 5 variables (Sand%, Clay%, OM%, Db, and ratio of CEC to clay), with 12 coefficients. The corresponding R2 values ranged from about 0.15 to 0.9, by soil group. The reported RMSE values by Børgesen and Schaap (2005) for PTF water-retention estimates at −10 and −100 kPa ranged from about ±3.3 to ±5.9%, by volume.

3.3.

Permanent wilting point

The permanent wilting point should generally increase with increasing field capacity, clay and OM content. The best-fitted linear predictor equation for the NB and NS values confirmed this, i.e.,

The scatter plot about the best-fitted regression line was found to be slightly curvilinear. In addition, this equation would not be consistent with the following requirements:

PWPWsoil = 0.23 × CLAYWsoil + 0.33 × OMWsoil + 0.18 × FCWsoil

• that FC values should always be positive, regardless of specific sand, OMWsoil and SPWsoil values; • that FCWsoil values should always be less than SP; • that FCWsoil should approach 0 as SP approaches 0; • that FCWsoil should decrease with increasing sand content.

As to be expected, PWPWsoil is strongly affected by increasing levels of internal soil surface area, as afforded by increasing organic matter and clay content. This equation, however, is not consistent with the expectation that all PWP values need to be
The nonlinear equation that contained SANDWsoil , OMWsoil and SPWsoil as predictor variables and ensured that 0 ≤ FCWsoil ≤ SPWsoil for all cases was formulated and fitted as follows: FCWsoil

(R2 = 0.692)

PWPWsoil

(8)



−0.511 × CLAY

= FCWsoil 1−exp





=SPWsoil 1−exp

(R2 = 0.958)



−0.588(1 − SANDWsoil ) − 1.73 × OMWsoil SPWsoil

Wsoil

− 0.865 × OMWsoil



FCWsoil

2

(R = 0.651)

(9)

For small CLAY and OM values, this equation reduces to (7)

With this equation, the best-fitted R2 value dropped slightly, but the resulting scatter plot of actual data versus best-fitted values became linear, as shown in Fig. 2, with the field capacity being quite low in sandy soils, as to be expected. Specifically, this equation reduces to FCWsoil ≈ 0.588(1 − SANDWsoil ) + 1.73OMWsoil when SANDWsoil approaches 1 and OMWsoil approaches 0, i.e., FC decreases with increasing sand content, and increases with increasing OM content. Shown in Fig. 2 (right) are 3D visualisations of how the best-fitted FCWsoil vary with changing SANDWsoil and OMWsoil in the NB and NS soils. With the above equation, FCWsoil can now be estimated within an MAE of ±0.032 g g−1 for the entire range of the varying soil conditions, or an RMSE of ±0.048 g g−1 . In comparison, Schaap et al. (1999) reported an RMSE of ±0.107 cm3 cm−3 for volume-based field capacity estimates using categorical texture specifications produced, while specifying the corresponding SAND, SILT, CLAY values produced an RMSE of

PWPWsoil ≈ 0.511 × CLAYWsoil + 0.865 × OMWsoil

(10)

i.e., this equation is consistent with the expectation that PWP increases with increasing CLAYWsoil and OMWsoil content. Shown in Fig. 3 (left) is the scatter plot of the source data pertaining to PWPWsoil versus their best-fitted values. This plot is again linear, and includes the “outlier” for peat. Also shown in Fig. 3 are the resulting 3D plots that show how PWPWsoil varies with clay, OMWsoil , and FCWsoil . In terms of precision, estimated PWPWsoil values were found to have an MAE of ±0.026, and an RMSE of ±0.035 across the varying soil conditions. In comparison, the above-mentioned Pachepsky and Rawls (1999) study reported an average RMSE of ±0.017 m3 m−3 and an RMSE range from about ±0.006 to ±0.024 m3 m−3 for soil moisture retention at −1500 kPa, using a compounding combination of linear and quadratic terms for the final suggested PTF. The corresponding R2 values ranged from about 0.63 to 0.98, by soil group. Børgesen and Schaap (2005) reported point and parametric PTF RMSE values water retention at −1500 kPa ranging from about ±0.023 to ±0.036 m3 m−3 .

307

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

Fig. 3 – Left: permanent wilting point (PWPWsoil ) by weight, in %, according to Eq. (9). Top: predicted PWPWsoil versus actual data (outlier on top right refers to peat). Center and bottom: 3D visualisations of how PWPWsoil varies with increasing FCWsoil , OM fraction and clay fraction. Right: saturated hydraulic conductivity (log10 Ksat ), in cm h−1 , according to Eq. (12). Top: predicted log10 (Ksat ) versus actual data. Center and bottom: 3D visualisations of how log10 (Ksat ) varies with increasing bulk density, sand fraction and OM fraction.

3.4.

Saturated hydraulic conductivity (Ksat , in cm h−1 )

Saturated hydraulic conductivities should increase with increasing sand content because sandy soils have large pores, and should decrease with increasing clay content. More generally, hydraulic conductivities should decrease logarithmically with decreasing pore space, or with increasing soil bulk density. The best-fitted linear predictor equation for the NB and NS values for log10 (Ksat ) confirmed these expectations, as follows: log10 (Ksat ) = 3.5 − 2.8 × Dbsoil + 2.1 × SANDWsoil

(R2 = 0.791) (11)

This equation, however, would be inconsistent with the expectation that Ksat should be 0 when Dpsoil = Dbsoil . Reflecting on the Ksat pore-space restrictions caused by Dpsoil − Dbsoil led to the following best-fitted formulation: log10 (Ksat soil ) = −0.98 + 7.94 log10 (Dpsoil − Dbsoil ) + 1.96 × SANDWsoil

(R2 = 0.795)

(12)

According to this, Ksat values are estimated to vary from 0 (as in solid rock, or fully compacted soils), to about 3710 cm day−1 (155 cm h−1 ) for a mineral soil (with Dpsoil = 2.65

308

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

Table 3 – Best-fitted results for Db (g cm−3 ), FCWsoil (g g−1 ), PWPWsoil (g g−1 ), and Ksat (cm h−1 ) including their respective a, b, c, and d coefficients, R2 , and RMSE for the New Brunswick and Nova Scotia soils data, using Eqs. (13)–(16) Property

ax

bx

cx

dx

R2

0.022 −0.004

6.12 −0.79

0.831

0.14

0.179

MAE

RMSE

Dbsoil

1.173 −0.049

0.83 −0.08

FCWsoil

0.588 ±0.016

1.734 ±0.049

0.958

0.032

0.048

PWPWsoil

0.511 ±0.025

0.865 ±0.057

0.651

0.026

0.035

Log10 (Ksat soil )

−0.98 ±0.11

7.94 ±0.48

0.795

0.381

0.488

1.96 ±0.21

Subscript x for a, b, c and d means Db, FC, PWP, or Ksat , as pertinent by row.

and Dbsoil = 1.23 g cm−3 at the soil surface, no sand), and to about 7.3 cm day−1 (0.3 cm h−1 ) for an organic soil (with Dpsoil = 1.3 and Dbsoil = 0.157 g cm−3 at the soil surface). Shown in Fig. 3 (right) are (i) the resulting scatter plot of fieldobserved versus best-fitted values, and (ii) 3D plots of how Ksat varies with increasing bulk density, and increasing OM and sand content. In all of this, the expected effect of OMWsoil on Ksat is indirect by way of Dbsoil and Eq. (5). The best-fitted RMSE value (±0.49) for log10 (Ksat ) is slightly better or similar to the corresponding RMSE values determined by Schaap et al. (1999): these authors used texture class (RMSE = ± 0.63), sand, silt, clay% (RMSE = ±0.60), and sand, silt, clay% plus soil bulk density (RMSE = ±0.533) to estimate log10 (Ksat ) by various means, including multiple regression and neural network analysis. Børgesen et al. (2007) reported best-fitted log10 (K) RMSE values (parametric PTFs, 18 formulations) from ±0.46 to ±0.78 for three Danish soil data set (n = 1618, 1609 and 152 soil layers, each with 4–7 water content readings from −1 to −1500 kPa), and K determinations from 0 to −80 h Pa). The Ksat numbers estimated for organic soils with Eq. (12) compare well with the field-estimated average by Schlotzhauer and Price (1999), namely 15 ± 8.1 RMSE cm day−1 , or 0.63 ± 0.34 RMSE cm h−1 . Furthermore, the best-fitted MAE and RMSE values in Table 3 suggest that actual Ksat values may be 2.4–3 times higher or lower than the estimated average. For organic soils, predicted Ksat values could then be as low as 2.4, or as high as 22 cm day−1 . Actual values for Ksat in well humified and undisturbed peat soils were in fact reported to

varied from 4.4 to 33 cm day−1 (Paivanen, 1973; Schlotzhauer and Price, 1999).

3.5.

Summary for the NB and NS equations

We now summarize and generalize the above expressions for the NB and NS soils as follows: Dbsoil =

aDb + (Dpsoil − aDb − bDb SANDWsoil )[1 − exp(−cDb DEPTH)] 1 + dDb OMWsoil (13)

FCWsoil



= SPWsoil

 1 − exp



−aFC (1 − SANDWsoil ) − bFC OMWsoil SPWsoil

(14) PWPWsoil



= FCWsoil 1 − exp

 −a

PWP CLAYWsoil

− bPWP OMWsoil



FCWsoil (15)

log10 (Ksat )=aKsat + bKsat log10 (Dpsoil − Dbsoil ) + cKsat SANDWsoil (16) Altogether, there are 4 adjustable coefficients for estimating Db, 2 for FC and PWP each, and 3 for Ksat . The best-fitted values for these coefficients are listed in Table 3, along with

Table 4 – Correlation coefficients between the a, b, c and d parameters of Eqs. (13)–(16), New Brunswick and Nova Scotia soil data Dbsoil ax bx cx dx

FCWsoil ax bx

ax

bx

1 −0.007 −0.45 0.77

1 0.75 0.062

cx

dx

1 −0.16

Log10 (Ksat soil )

ax

bx

cx

ax bx cx

1 −0.043 −0.89

1 −0.15

1

1

ax

bx

PWPWsoil

ax

bx

1 −0.58

1

ax bx

1 −0.34

1

309

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

Table 5 – FC and PWP model comparisons, using NB, NS and ISRIC (spodosol and histisol) data a

b

FC NB and NS soil data ISRIC spodosol and histisol data ISRIC spodosol and histisol data, optimized NB, NS + ISRIC (spodosol and histisol) data NB, NS + ISRIC (spodosol and histisol) data optimized

0.59 0.59 0.59 0.59 0.50 ± 0.03

1.74 1.74 1.74 1.74 2.73 ± 0.25

PWP NB and NS soil data ISRIC spodosol and histisol data NB, NS + ISRIC (spodosol and histisol) data

0.511 0.511 0.511

0.865 0.865 0.865

their approximated standard errors of estimate. Also in this table are the values for R2 , MAE, and RMSE, for each of the four regression models. The extent of inter-parametric correlation among the various a, b, c and d coefficients is shown in Table 4. These numbers should be as close to zero as possible to narrow the equifinal solution space for these best-fitted coefficients. For example, the −0.89 correlation between aKsat and cKsat implies that an increase in cKsat will produce a corresponding decrease in aKsat . Hence, larger non-zero inter-coefficient correlations numbers imply larger uncertainties about the best-fitted coefficient values. For Eqs. (13) to (16), the achievable precision varied in the following order: FC > PWP > Db > Ksat . This order is likely influenced by the extent to which soil structure (or state of soil aggregation) affects these variables. In general, soil structure should have (i) some effects on the FC and PWP pressure-plate determinations, (ii) substantial effects on in situ Db values as affected by soil depth and OMWsoil , and undisturbed aggregate size in general, and (iii) large effects on Ksat , on account of disproportionate rates of flow through pores, root channels and cracks of varying size. With organic soils, Ksat is also strongly affected by the varying degrees of humification, with well-humified matter being much less permeable than fibrous matter (Pepin et al., 1992; Paquet et al., 1993). To account for this with Eq. (12) and for organic soils only, we set aKsat = (2.5 ± 0.2) − (0.046 ± 0.004) vanPost index2 (R2 = 0.89), based on the peat humification data compiled by Letts et al. (2000). From this, it follows that the aKsat value of −0.97 in Eq. (12) would correspond to a van Post humification index of about 8.5. We further note that the aKsat value of 2.5 ± 0.2 would then not only correspond to a van Post humification index of 0, but would also reflect reported flow rates through saturated and freshly milled sphagnum peat (Shibchurn et al., 2005).

n

R2

MAE

RMSE

1 1 1.064 ± 0.008 1.064 1.058 ± 0.009

215 233 233 448 448

0.921 0.650 0.797 0.857 0.877

0.031 0.061 0.057 0.0443 0.0406

0.044 0.095 0.075 0.0625 0.0584

1 1 1

190 236 426

0.65 0.58 0.60

0.026 0.043 0.036

0.035 0.058 0.048

f

applying Eqs. (14) and (15) – are summarized in Table 5. In combining the two data sets, we addressed the procedural difference between the NB, NS and ISRIC data for Db (converted to SP): the latter refer to soil cores equilibrated at a soil moisture potential of −33 kPa, while the former refer to undisturbed soil cores. In Eq. (14), this difference was accounted for by replacing f SPWsoil with SPWsoil , with f as an adjustable exponent, which then led to the best-fitted value of f = 1.064 ± 0.008. As indicated in Table 5, R2 values between actual versus best-fitted FC values for the ISRIC data would have dropped from 0.797 to 0.650, if the f exponent would remain at its default value of 1. Adjusting the “a”, “b” and “f” coefficients for the combined NB, NS, and ISRIC data (f adjusted for the ISRIC data only) improved the best-fitted R2 , MAE and RMSE values, but increased the “b” coefficient from 1.74 to 2.73 while the “a” and “f” values nearly remained the same. Hence, there is a greater amount of uncertainty associated with the OM coefficient than with the sand and “f” coefficients. For PWP, Eq. (15) was used without modification for the ISRIC, NB and NS data, and the best-fitted results are also listed in Table 5. Actual versus best-fitted FC and PWP values for the combined NB, NS and ISRIC data (Table 5) are plotted in Fig. 4 (top and middle). These plots confirm that Eqs. (14) and (15) represent the FC and PWP variations of both data sets quite well.

3.7. Verifying the Ksat results for NB and NS soils with data from similar soils elsewhere Combining the Ksat -relevant data of the Universal Soil Database (UNSODA) with the NB and NS data as part of verifying Eq. (16) yielded: log10 (Ksat ) = (−1.05 ± 0.08) + (6.1 ± 0.4) log10 (Dpsoil − Dbsoil ) + (2.2 ± 0.1)SANDWsoil

3.6. Verifying the FC and PWP results for NB and NS with data from similar soils elsewhere

× (n=481; R2 =0.52; RMSE = 0.68; MAE = 0.54) (17)

Since the dominant NB and NS soil types belonged to spodosols and histisols, we selected the corresponding Db, SP, FC, PWP data for spodosols and histisols from the ISRIC database. We then analyzed the pool of combined data in three ways: (a) NB and NS soils alone, (b) ISRIC spodosol and histisol data alone, and (c) in combination. The resulting best-fitted equation coefficients, R2 , MAE and RMSE values – generated by

This shows that the above log10 (Ksat ) formulation remain valid in its general form, but the values of its three coefficients changed slightly. In this instance, the largest change was associated with the log10 (Dpsoil − Dbsoil ) coefficient, namely a drop from 7.9 (Table 3) to 6.1. This change may again relate to procedural data differences, e.g., using estimated Dpsoil values from

310

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

making a few further modifications to the above equations to account for the wider range of soil conditions improved the overall fit between actual and resulting best-fitted estimates, and generated the following expressions: Db =

aDb + [Dp −aDb −bDb (1 − CLAYWsoil )][1−exp(−cDb DEPTH)] 1+dDb OMWsoil +eDb Traffic;

(18)

] FCWsoil = SPWsoil [cFC + (dFC − cFC )CLAY0.5 Wsoil

 a SAND  FC Wsoil − bFC OMWsoil

exp −

(19)

SPWsoil



PWPWsoil = FCWsoil cPWP + (dPWP − cPWP )CLAY0.5 Wsoil



 a SAND  PWP Wsoil − bPWP OMWsoil

exp −

Fig. 4 – Scatterplots of ISRIC data for FCWsoil (top) and PWPWsoil (middle) versus the corresponding estimates from the NB- and NS-derived Eqs. (7) and (9) (ISRIC spodosols and histisols only, n = 234), in % by weight. Bottom: scatterplot of the UNSDODA data (n = 459) for log10 Ksat (bottom) versus the corresponding estimates from the NBand NS-derived Eq. (12) (n = 370). Plots also show actual data versus best-fitted NB and NS values, for comparison; Ksat units: cm h−1 .

known SPWsoil and Dbsoil values (NB and NS data) versus direct Dpsoil measurements (ISRIC data). The plot of actual versus best-fitted log10 (Ksat ) values in Fig. 4 (bottom) demonstrates the general conformance of the above log10 (Ksat ) formulation to the data from both sources.

3.8. Extending the generalized Db, FC, and PWP equations to other soils The ISRIC data were then used to extend the above formulation for Db, FC and PWP to soils other than histisols and spodosols. Altogether, 13,103 data rows of the ISRIC database (each row representing a particular soil layer) were used for this task. Starting with Eqs. (13)–(15), and

FCWsoil

(20)

Soils showing unusually high Db values were flagged by introducing the “traffic” variable, and by setting this variable equal to 1 when deemed appropriate, and 0 otherwise. The resulting best-fitted values for R2 , MAE, and RMSE and the a, b, c and d coefficients for each of these equations are listed in Table 6 for Db, and in Table 7 for FC and PWP. For the most part, the resulting R2 values varied widely, but the results for spodosols and histisols combined with temperate to boreal forest soils and wetlands were again fairly consistent with the results already listed in Tables 3 and 5. In general, the Db associated a, b, c, d and e coefficients were found to change by soil order. These changes appear to be systematic, and likely reflect direct influences of climate, vegetation, soil substrate, extent of soil weathering and land use on soil bulk densities. For example, vertisols tend to have very high bulk densities at the surface, especially when organic matter levels are low, likely due to extensive soil crusting. In contrast, forest soils with organic forest floor layers on top tend to have fairly low densities at the mineral soil surface. Andisols as a group show the greatest Db differences. Some of their Db values are similar to aridisols, entisols, and mollisols, but other Db values are unusually low throughout the profile, have a low dependence on organic matter content, and have no discernable dependence on clay content. Some of these differences are likely related to a limited degree of soil weathering of the volcanically derived soil substrates. In contrast, the very well weathered oxisols and ultisols show a fairly high bulk density with little overall dependence on soil depth. The strongest effect of clay on soil bulk density occurs with clay-enriched mollisols (e.g., argiaquolls, argiudolls), and with soils having a low Ca content (e.g., natrustolls). This is undoubtedly reflecting soil-related differences by clay type. We note that Eqs. (19) and (20) can be simplified in various ways. For example, the condition that SAND = CLAY = OM = 0 (or SILT = 1) suggests that FCWsoil = cFC SPWsoil

and

PWP = cPWP FCWsoil

311

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

Table 6 – Best-fitted results for soil bulk density (Db, g cm−3 , Eq. (18)), including optimal a, b, c, d, e parameters, R2 , RMSE, sample size (n), and relative error %, by soil order, with exceptions; ISRIC data RMSE %

R2

n

r

6.27 6.27

16.3 18.8

0.93 0.67

15 408

0.96 0.82

9.3 −1.9

0.022 0.010

6.27 6.27

18.6 19.4

0.41 0.15

1437 163

0.64 0.39

−0.4 7.0

1.10 1.10

0.022 0.002

6.27 6.27

13.8 9.6

0.19 0.72

990 14

0.44 0.85

−2.0 0.7

1.50

1.10

0.022

6.27

12.5

0.11

1007

0.33

0.0

General Exceptionsc

1.50 1.50

1.10 1.40

0.022 0.022

6.27 6.27

14.4 16.5

0.23 0.21

3487 802

0.48 0.46

−2.8 10.6

Alfisols

General Exceptionsd

1.65 1.18

1.10 0.89

0.022 0.022

6.27 6.27

13.8 20.1

0.11 0.29

3116 394

0.33 0.54

−3.8 −7.1

Oxisols

General Exceptionse

1.40 1.40

1.20 1.20

0.002 0.001

6.27 6.27

15.8 16.4

0.07 0.14

245 21

0.26 0.37

−0.4 −3.7

Ultisols

General Exceptionsf

1.55 1.55

1.00 1.00

0.005 0.001

6.27 6.27

13.4 16.6

0.11 0.16

1327 98

0.33 0.40

−1.6 −6.5

Vertisols

General

1.90

1.20

0.010

6.27

10.6

0.14

588

0.37

−3.5

Andisols

General Exceptionsg Exceptionsh

1.40 1.00 1.00

1.10 0.00 0.00

0.050 0.001 0.001

6.27 2 2

27.5 21.6 16.5

0.48 0.06 0.05

161 229 38

0.69 0.24 0.22

−6.4 −11.4 0.2

Soil order

Data groups

Surface a

Clay b

Histosols Spodosols

General General

1.18 1.18

0.89 0.70

0.022 0.022

Inceptisols

General Exceptionsa

1.18 1.50

0.89 1.10

Entisols

General Exceptionsb

1.50 1.50

Aridisols

General

Mollisols

a b c d e f g

h i

Depth c

OM d

Taffic e

0.2

0.5

Relative error %i

Dystropept, eutropept, halaquept, tropept, ustropept and xerumbrept. Cryorthent. Aquoll, argiaquoll, argiudoll, cryaquoll, duraquoll, natrustoll and paleustoll. Cryoboralf, fragiboralf, kandiudalf, kanhaplustalf and paleudalf. Kandiudox. Haploxerult, kandihumult, kanhaplustult, palehumult, rhodudult and udult. Haplaquands, hapludands, haplustands, melanudands, udivitrands, vitraquands, vitricryands, vitritorrands, vitrixerands, with measured Db < 1.4. Melanudands, udivitrands, vitraquands, vitricryands, vitritorrands, vitrixerands, with measured Db > 1.4. Relative error = mean error of predicted Db values/mean of actual Db values, by soil order.

Similarly, for a mixture of silt and organic matter, it follows that

Similarly, for pure clay, it follows that FCWsoil = dFC SPWsoil

PWP = dPWP FCWsoil

and

 b OM  FC Wsoil

FCWsoil = cFC SPWsoil exp −

For a mixture of silt and sand, one obtains

 a SAND  FC Wsoil

FCWsoil = cFC SPWsoil exp −

SPWsoil

and

 b  PWP OMWsoil

PWPWsoil = cPWP FCWsoil exp −

and

 a  PWP SANDWsoil

PWPWsoil = cPWP FCWsoil exp −

SPWsoil

FCWsoil

Simplifying the above FC expression to FCWsoil = cFC SPWsoil and FCWsoil = [cFC + (dFC − cFC )CLAY0.5 ]SPWsoil , and using these Wsoil

FCWsoil

Table 7 – Best-fitted results for FCWsoil and PWPWsoil , including regression coefficients, R2 , MAE, RMSE, and sample size (n); ISRIC data Property FCWsoil

n

bx

13,103 13,103 13,103

PWPWsoil

ax

0.103 ± 0.057

0.785 ± 0.009

13,088 13,088 13,088

0 ± 0.52

−1.4 ± 8.3

cx

dx

R2

MAE

RMSE

0.678 ± 0.002 0.432 ± 0.003

0.963 ± 0.004

0.717 0.814

4.95 4.01

6.86 5.56

0.565 ± 0.011

0.991 ± 0.008

0.832

3.96

5.49

0.518 ± 0.002 0.132 ± 0.015

0.818 ± 0.006

0.518 0.725

3.95 2.88

5.37 3.91

0.170 ± 0.010

0.832 ± 0.003

0.725

2.84

3.86

Equation FCWsoil = cx SPWsoil FCWsoil = [cx + (dx − cx ) ]SPWsoil CLAY0.5 Wsoil Eq. (19) PWPWsoil = cx FCWsoil PWPWsoil = [cx+ (dx − cx ) CLAY0.5 ]FCWsoil Wsoil Eq. (20)

312

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

simplifications to model all of the compiled ISRIC data, generated R2 values of 0.717 and 0.814, respectively. Similarly, simplifying PWP to PWPWsoil = cPWP FCWsoil and PWPWsoil = [cPWP cFC + (dPWP − cPWP )CLAY0.5 ]FCWsoil generated R2 values Wsoil of 0.518 and 0.725, respectively (Table 7). All of this implies that • pore space alone is the most dominant factor in determining moisture retention at field capacity; • specifying clay content adds significantly to explaining the observed FC variations (Hillel, 1982; Yong and Warkentin, 1975); • organic matter effects on FC are mainly related to the influence of organic matter on soil bulk density, or pore space; apart from this, organic matter may also have additional effects on FC, presumably on account of organicmatter induced water repellency, especially under dry and arid conditions, or the lack thereof under humid and wet soil conditions (Ma’shum and Farmer, 1985; Neinhuis and Barthlott, 1997); • most of the ISRIC PWP data variations can be attributed to changes in FC alone, with clay once again adding significantly to the overall observed PWP variations; in contrast, effects of sand and organic matter do not appear to influence PWP significantly, other than through their influence on soil bulk density. Applying Eqs. (19) and (20) with the best-fitted ISRICderived coefficient values in Table 7 to the NB and NS data produced the following results • for FCWsoil : R2 = 0.815, MAE = 2.1%, RMSE = 3.11 (by weight), and using a best-fitted exponent value of f = 0.901 ± 0.004 f after substituting SPWsoil with SPWsoil ; • for PWPWsoil : R2 = 0.39, MAE = 5.6, RMSE = 6.6 (by weight). These results are further illustrated in Fig. 5, which shows a visual comparison of the actual versus best-fitted FC and PWP

Fig. 5 – Scatterplots for actual and best-fitted ISRIC FCWsoil and PWPWsoil data values (%, n = 13,103; Eqs. (19) and (20)), with actual data versus best-fitted NB and NS values overlain.

values for the NB and NS data in reference to the corresponding ISRIC FC and PWP scatter plots. This comparison reveals that the best-fitted results generated from the more universally derived Eqs. (19) and (20) conform to the results derived with Eqs. (14) and (15), but are not as accurate. Hence, calibrat-

Fig. 6 – Variation of soil depth, pore space (SPWsoil ), soil solid space, FCWsoil , PWPWsoil , potentially available water and Ksat calculated for a soil with increasing clay fraction, containing complementary portions of sand, with a silt fraction set at 0.25 by weight, and organic matter fractions set at 0 (left) and 0.1: under natural conditions (left, middle) and recompressed to the same density of the purely mineral soil (right). Based on Eqs. (2), (7)–(9) and (12).

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

ing the equations over a wider range of soil conditions likely masks some of the details that affect soil water retention, soil density and soil permeability regionally and locally. Across varying regions and soil climates, the above ISRIC-derived FC and PWP equations appear to be fairly robust, especially in terms of their simplified expressions, thereby indicating that a considerable portion of soil moisture retention can be estimated from simple pore space and clay specifications alone.

3.9.

Application

The above results can now be used to estimate effects of organic matter additions and soil compaction on Db, FC, PWP and Ksat for any particular texture condition. Doing this with the Table 3 entries and Eqs. (2), (7)–(9) and (12) is illustrated in Fig. 6 for a soil with increasing clay content (or decreasing sand content) and a fixed silt content (25% by weight) with 0 and 10% organic matter content (by weight) under natural conditions, and when the soil with the 10% organic matter content is re-compressed. These results show that • an increased clay content would decrease Ksat strongly while FCWsoil and PWPWsoil would both increase at the same time; the amount of maximum available water (=FCWsoil − PWPWsoil ), however, would decrease; • adding organic matter from 0 to 10% (by weight) would double the soil volume, and would increase the potentially available water for the same soil even more, at any clay content; this would occur in spite of a small but persistent decrease in overall soil depth (or increase in soil bulk density) with increasing clay content at 10% organic matter; whether the slight decrease in soil depth with increasing clay content but constant organic matter content at 10% is realistic remains to be seen, and would likely be affected by the extent and timing of the aggregation processes within the soil; • re-compacting an OM-amended soil to its original density would decrease Ksat below its non-amended level; on compaction, the added OM would be squeezed into the pore space that the non-amended soil would have had, and water would then flow considerably more slowly through the recompacted soil than the non-amended soil.

4.

Concluding remarks

To summarize, hydrological modelling uncertainties following the above PTF calibrations generally fall into an RSME range of ±3 to ±6% for pressure-plate determined SP, FC, and PWP values, increase to about ±14% for in situ Db values, i.e., natural soil conditions for which Db generally increases with increasing soil depth, and remain fairly large for log10 (Ksat ) with best-fitted RSME values from about ±0.4 to ±0.6. As intended, this approach produces physically feasible and generalizable estimates for Db, FC, PWP and Ksat from soil texture, organic matter content and soil depth specifications, as these would vary from highly porous organic soils to highly compacted clay soils. In addition, the approach reveals how the best-fitted regression coefficients for the ISRIC Db data tend to vary systematically by soil order or soil order groups. In terms

313

of application, the above equations, together with a similar derivations for soil thermal conductivity and soil mechanical resistance, can be used for initializing hydrological models that are designed to simulate effects of weather, climate and land-use changes on hydro-thermal, -mechanical and biological soil conditions as these vary daily and year-round (see, e.g., Balland and Arp, 2005; Balland et al., 2006; Vega-Nieva et al., 2009). To this effect, special considerations need to be given to track local changes in soil bulk density for proper model initialization, because these may vary with changing surface traffic, weather, and soil depth. Weather-induced changes in soil density and soil permeability refer to shrinking and swelling, and freezing and thawing. Once these changes are known, the above derived expressions can be expected to provide reasonable estimates for SP, FC, PWP and Ksat for a wide range of soil conditions in a systematic way. Local to regional calibrations should further enhance the application potential and precision of the above approach, from saturated to unsaturated soil conditions.

Acknowledgements This research was supported by the NSERC Discovery Grants to PAA, the Nexfor/Bowater Forest Watershed Research Centre at UNB, by Environment Canada in support of the continuing development of the Forest Hydrology Model ForHyM (c/o T. Clair), and the Georgia Basin Project. This work was further supported by a soil trafficability research contract with DND, Canada, and by NSERC’s Sustainable Forest Management Network on modeling and mapping hydrologically sensitive areas. We give special thanks to the reviewers of this manuscript, and to Jose Miguel Soria Ugalde from the University of Guanajuato (Mexico) for his review of an earlier version of this manuscript, and to Lee Lihui and Cheng-fu Zhang for manuscript assistance.

Appendix A. Definitions and basic equations There are three key soil moisture levels (Hillel, 1982): SP: the soil is at saturation when all the pores are filled with water. FC: after water has drained from the soil due to gravity, the larger pores are drained, but the smaller pores remain filled with water; at this point, the soil is said to be at field capacity; generally, this occurs at a soil matric potential of −33 kPa. PWP: plants remove water from the soil until they begin to wilt, but even at the permanent wilting point substantial amounts of soil water remain, especially in organic and clayey soils. However, this water is held too tightly to permit absorption by plant roots. Generally, this occurs at a soil matric potential of −1500 kPa. In reference to Fig. 1, the various terms of the above equations are defined as follows: SILTWsoil =

Wsilt Wparticles

(A.1)

314

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

SANDWsoil =

CLAYWsoil =

Wsand Wparticles Wclay Wparticles

(A.2)

(A.4)

Including coarse fragments as part of the whole requires a re-calculation of all of the above soil volume/weight relationships. Formally, CF volume and weight fractions are defined by CFV = VCF /Vtotal and CFW = WCF /Wsollid. The density of CF is given by DCF = WCF /WCF . The bulk density of the soil that includes CF, i.e., Db, can then be obtained from

(A.5)

1 − CFW CFW 1 + = Db DCF Dbsoil

(A.3)

MINERALSWsoil = Wminerals /Wparticles = SANDWsoil + SILTWsoil + CLAYWsoil

OMWsoil =

WOM Wparticles

A.1. Dealing with coarse fragments

(A.12)

if CFW is known. Alternatively, if CFV is known, Two of the commonly used soil density measures are particle density (Dp) and bulk density (Db). For the purpose of this paper, and in keeping with common practice, these densities do not include moisture. With CF included, we use: Wsolid Vsolid

(A.6)

W Db = solid Vtotal

(A.7)

Dp =

Db = Dbsoil (1 − CFV ) + CFV DCF

(A.13)

The average particle density of the soil plus CF can be calculated from CFW 1 − CFW 1 = + Dp DCF Dpsoil

(A.14)

From this, the total volumetric saturation point SPV can then be obtained from

Without CF, we set Dpsoil =

Dbsoil =

Wparticles Vparticles Wparticles Vsoil

SPV = (A.8)

(A.9)

Soil bulk density, particle density and pore space are all related to the above quantities, i.e., Vpores /Vtotal = 1 − Db/Dp. Similarly, Vpores /Vtotal = 1 − Dbsoil /Dpsoi . In terms of OM and soil minerals, we set DOM = WOM /VOM = 1.3 g cm−3 for most organic materials, and Dminerals = Wminerals /Vminerals = 2.65 g cm−3 for most silicate and carbonate minerals. From these particular definitions, it can be shown that OMWsoil 1 − OMWsoil 1 = + Dpsoil DDOM DMIN

(A.10)

Soil water contents in term of weight or volume fraction (subscript “V” for volume) are given by WATERvsoil = Vwater /Vsoil WATERWsoil = Wwater /Wparticles , and WATERV = Vwater /Vtotal . Converting soil moisture volume to weight fraction is readily achieved by setting WATERWsoil = WATERVsoil (Dwater /Dbsoil ), where Dwater is the density of water, equal to 1 g cm−3 . We now define SP, FC, and PWP as weight and volume fractions using Eqs. (4)–(6), with the following symbols: SPWsoil , FCWsoil , PWPWsoil , SPVsoil , FCVsoil , PWPVsoil , SPV , FCV , PWP. The rate of water flow through soils at saturation is determined by the hydraulic conductivity K, defined by Darcy’s law: K=

A(H/L) dQ/dt

(A.11)

where H is the total hydraulic head difference from the top to the bottom of a saturated soil column, L is the length of this column, dQ/dt is the rate of water flowing through this column, and A is the cross-sectional area.

Vpores Db =1− Vtotal Dp

(A.15)

Since SPV , FCV and PWPV are all reduced by the same proportion with increasing CF, it follows that the ratio of volumetric field capacity over volumetric pore space must remain the same regardless of CF. Hence, and FCV = (FCVsoil /SPVsoil )SPV PWPV = (PWPVsoil /SPVsoil )SPV . With respect to soil percolation, all of the above remains valid for the soil that fills the space among the coarse fragments, but with Darcy’s law accounting for the presence of coarse fragments through a proportional adjustment of the area through which the flow occurs, i.e.,

 dQ  dt

=K

 H  L

(A − ACF )]

(A.16)

where ACFb is the effective cross-sectional area occupied by CF.

references

Agnew, J.M., Leonard, J.J., Feddes, J., Feng, Y., 2003. A modified air pycnometer for compost air volume and density determination. Can. Biosyst. Eng. 45, 627– 635. Balland, V., Arp, P.A., 2005. Modeling soil thermal conductivities over a wide range of conditions. J. Eng. Environ. Sci. 4, 549–558. Balland, V., Bhatti, J.S., Errington, R., Castonguay, M., Arp, P.A., 2006. Modeling soil temperature and moisture regimes in a jack pine, black spruce and aspen forest stand in central Saskatchewan (BOREAS). Can. J. Soil Sci. 86, 203–217. Blake, G.R., Hartge, K.H., 1986. Bulk density. In: Klute, A. (Ed.), Methods of Soil Analysis. Part 1. Physical and Mineralogical Methods, 2nd ed. Agron. Monogr. 9. ASA and SSSA, Madison, WI, pp. 363–382. Bouma, J., 1989. Using soil survey data for quantitative land evaluation. Adv. Soil Sci. 9, 177–213.

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

Børgesen, C.D., Schaap, M.G., 2005. Point and parameter pedotransfer functions for water retention predictions for Danish soils. Geoderma 127, 154–167. Børgesen, C.D., Iversen, B.V., Jacobsen, O.H., Schaap, M.G., 2007. Pedotransfer functions estimating soil hydraulic properties using different soil parameters. Hydrol. Process. 251, 123–150. CANSIS, 2000. Canadian Soil Information System. Available from http://www.sis.agr.gc.ca/cansis/. Fahmy, S.H., Rees, H.W., 1989. Soils of the Woodstock-Florenceville area, Carleton County, New Brunswick. Report No. 14. Edited by Research Branch, Agriculture Canada. pp. 33–51. Futter, M.N., Butterfield, D., Cosby, B.J., Dillon, P.J., Wade, A.J., Whitehead, P.G., 2007. Modeling the mechanisms that control in-stream dissolved organic carbon dynamics in upland and forested catchments. Water Resource Res. 43, W02424, doi:10.1029/2006WR004960. Gijsman, A., Jagtap, S.S., Jones, J.W., 2003. Wading through a swamp of complete confusion: how to choose a method for estimating soil water retention parameters for crop models. Eur. J. Agron. 18, 75–105. Heuscher, S.A., Brandt, C.C., Jardine, P.M., 2005. Using soil physical and chemical properties to estimate bulk density. Soil Sci. Soc. Am. J. 69, 51–56. Hillel, D., 1982. Introduction to Soils Physics. Academic Press, Inc./Harcourt Brace Jovanovich Publishers, pp. 155–166. Langmaid, K.K., MacMillan, J.K, Losier, J.G., 1976. Soils of Northern Victoria County, New Brunswick. Report No. 7. Edited by Research Branch, Canada Department of Agriculture. pp. 220–255. Leij, F.J., Alves, W.J., van Genuchten, M.Th., Williams, J.R., 1996. The UNSODA Unsaturated Soil Hydraulic Database; User’s Manual, Version 1.0. EPA/600/R-96/095, National Risk Management Laboratory, Office of Research and Development, U.S. Environmental Protection Agency, Cincinnati, OH. http://www.ussl.ars.usda.gov/models/unsoda.HTM. Letts, M.G., Roulet, N.T., Comer, N.T., Skarupa, M.R., Verseghy, D.L., 2000. Parametrization of peatland hydraulic properties for the Canadian land surface scheme. Atmosphere-Ocean 38, 141–160, doi:0705-5900/2000/0000-0141. Ma’shum, M., Farmer, V.C., 1985. Extraction and characterisation of water-repellent materials from Australian soils. Aust. J. Soil Res. 23, 623–626. Mellerowicz, K.T., Rees, H.W., Chow, T.L., and Ghanem, I., 1993. Soils of the Black Brook Watershed, St-André Parish, Madawaska County, New Brunswick. Edited by Research Branch, Agriculture Canada. pp. 26–37. Minasny, B., Hopmans, J.W., Harter, T., Eching, S.O., Tuli, A., Denton, M.A., 2004. Neural networks prediction of soil hydraulic functions for alluvial soils using multistep outflow data. Soil Sci. Soc. Am. J. 68, 417–429. Neinhuis, C.W., Barthlott, 1997. Characterisation and distribution of water-repellent, self-cleaning plant surfaces. Ann. Bot. 79, 667–677. Nemes, A., Schaap, M., Wösten, H., 2003. Functional evaluation of pedotransfer functions derived from different scales of data collection. Soil Sci. Soc. Am. J. 67, 1093–1102. Pachepsky, Y.A., Timlin, D., Varallyay, G., 1996. Artificial neural networks to estimate soil water retention from easily measurable data. Soil Sci. Soc. Am. J. 60, 727–733. Pachepsky, Y.A., Rawls, W.J., 1999. Accuracy and reliability of pedotransfer functioned as affected by grouping soils. Soil Sci. Soc. Am. J. 63, 1748–1757. Pachepsky, Y.A., Rawls, W.J., Gimenez, D., 2001. Comparison of soil water retention at field and laboratory scales. Soil Sci. Soc. Am. J. 65, 460–462. Pachepsky, Y.A., Rawls, W.J., 2004. Development of pedotransfer functions in soil hydrology. Dev. Soil Sci., 30.

315

Paivanen, J., 1973. Hydraulic conductivity and water retention in peat soils. Report from the Faculty of Agriculture and Forestry of the University of Helsinki, Finland. Acta Forestalia Fennica 129, 1–70. Paquet, J.M., Caron, J., Banton, O., 1993. In-situ determination of the water desorption characteristics of peat substrates. Can. J. Soil Sci. 73, 329–339. Pepin, S., Plamondon, A., Stein, J., 1992. Peat water content measurement using time domain reflectometry. Can. J. Forest Res. 22, 534–540. Rawls, W.J., Brakensiek, D.L., 1989. Estimation of soil retention and hydraulic properties. In: Morel-Seytoux (Ed.), Unsaturated Flow in Hydrologic Modeling, pp. 275–300. Rawls, W.J., Gish, T.J., Brakensiek, D.L., 1991. Estimation of soil water retention from soil properties—a review. Adv. Soil Sci. 16, 213–235. Rawls, W.J., Pachepsky, Y.A., Ritchie, J.C., Sobecki, T.M., Bloodworth, H., 2003. Effect of soil organic carbon on soil water retention. Geoderma 116, 61–76. Rees, H.W., Langmaid, K.K., Losier, J.G., Veer, C., Wang, C., Wells, R.E., Fahmy, S.H., 1992. Soils of the Chipman-Minto-Harcourt region of New Brunswick. Report No. 11. Edited by Research Branch, Agriculture Canada. pp. 209–314. Ritchie, J.T., Ratliff, L.F., Cassel, D.K., 1987. Soil Laboratory Data, Field Descriptions and Field Measured Soil Water Limits for some Soils of the United States. ARS Technical Bulletin, ARS Agric. Research Service and Soil Conservation Service/State Agric. Exp. Stations, Temple, TX. Saxton, K.E., Rawls, W.J., Romberger, J.S, Papendick, R.I., 1986. Estimating generalised soil-water characteristics from texture. Soil Sci. Soc. Am. J. 50, 1031–1036. Schaap, M.G., Bouten, W., 1996. Modeling water retention curves of sandy soils using neural networks. Water Resource Res. 32, 3033–3040. Schaap, M.G., Leij, F.J., 1998. Using neural networks to predict soil water retention and soil hydraulic conductivity. Soil Tillage Res. 47, 37–42. Schaap, M.G., Leij, F.J., van Genuchten, 1999. Development of Pedotransfer Functions and Related Computer Programs. U.S. Salinity Laboratory, Riverside, CA, p. 7. Schaap, M.G., Nemes, A., van Genuchten, 2004. Comparison of models for indirect estimation of water retention and available water in surface soils. Vadose Zone J. 3, 1455–1463. Schaap, M.G., 2006. Models for indirect estimation of soil hydraulic properties. Encyclopedia Hydrol. Sci., doi:10.1002/0470848944.hsa078. Schlotzhauer, S.M., Price, J.S., 1999. Soil water flow dynamics in a managed cutover peat field: field and laboratory investigations. Water Resources Res. 12, 3675–3683. Shibchurn, A., Van Geel, P.J., Kennedy, P.L., 2005. Impact of density on the hydraulic properties of peat and the time domain reflectometry (TDR) moisture calibration curve. Can. Geotech. J. 42, 279–286. Soil Survey Laboratory Methods Manual, 2004. Soil Survey Laboratory Investigations Report No. 42. http://soils.usda.gov/technical/lmm/. Tamari, S., Wösten, J.H.M., Ruiz-Suarez, J.C., 1996. Testing an artificial neural network for predicting soil hydraulic conductivity. Soil Sci. Soc. Am. J. 60, 1732–1741. Tempel, P., Batjes, N.H., van Engelen, V.W.P., 1996. IGBP-DIS Soil Data Set for Pedotranfer Function Development. International Soil Reference and Information Center. http://www.isric.org/NR/exeres/545B0669-6743-402B-B79ADBF57E9FA67F.htm. Vega-Nieva, D.J., Castonguay, M., Ogilvie, J., Arp, P.A., 2009. Development of a modular terrain model to estimate daily variations in machine-specific forest soil trafficability, year-round. Can. J. Soil Sci., in press.

316

e c o l o g i c a l m o d e l l i n g 2 1 9 ( 2 0 0 8 ) 300–316

Webb, K.T., Thompson, R.L., Beke, G.J., Nowland, J.L., 1991. Soils of Colchester County, Nova Scotia. Report No. 19. Edited by Research Branch, Agriculture Canada. pp. 143–189. Wösten, J.H.M., Pachepsky, Y.A., Rawls, W.J., 2001. Pedotransfer functions: bridging the gap between available basic soil data and 2.0. J. Hydrol. 251, 151–162. Yong, R.N., Warkentin, B.P., 1975. Soil Properties and Behaviour. Elsevier Scientific Publ. Co., Amsterdam.

Yu, Z.C., Vitt, D.H., Campbell, I.D., Apps, M.J., 2003. Understanding Holocene peat accumulation pattern of continental fens in western Canada. Can. J. Bot. 81, 267–282. Zhang, C.F., Meng, F.-R., Bhatti, J.S., Trofymow, J.A., Arp, P.A., 2008. Modeling forest leaf-litter decomposition and N mineralization in litterbags, placed across Canada: A 5-model comparison. Ecol. Model. 219, 342–360.