Please see PDF version

(Wilconxon Prize 1989)

ACE Guided Transformation Method for Estimation of the Coefficient of Soil Water Diffusivity

Richard D. De Veaux and J. Michael Steele

Program in Statistics and Operations Research
Princeton University
Princeton, NJ 08544

Dataanalytic tools for choosing transformations to increase linear association are applied to a basic problem of soil physics, the determination of the coefficient of soil-water diffusivity. Data on Manawatu sandy loam illustrate the decisions the analyst must face and the quality of the estimates that the analyst can expect.

KEY WORDS: Bulge rule; Diffusion equation; Estimation of derivatives and integrals; Transformation choice.


1 1

i i

In many problems of science and engineering, one needs to estimate derivatives and integrals of functions that are determined empirically. One problem that requires such estimates and that seems particularly interesting from a statistical perspective is the estimation of the coefficient of soilwater diffusivity D(O). Despite being relatively unfamiliar to statisticians, this problem is a basic one of soil physics, and empirical findings are regularly reported. The feature of the problem that makes it unusual is that estimation of D(O) typically requires determining both derivatives and integrals of a function that is indirectly observed.
The present approach to the estimation of D(O) is guided by the view that statistical methods for dealing with data that exhibit stronglinear associations are well developed; consequently, many nonstandard problems are best addressed by transforming the data to achieve increased linear association. The fundamental notion of straightening curves and plots is well known in the statistical literature. Still, when that simple notion is applied to a technological problem that has traditionally been addressed by ad hoc means, the outcome can be surprising. In this instance, one finds a natural estimation process that has several advantages over methods that are widely applied. The resulting approach also has the interesting feature of incorporating dataanalytic exploration and conventional statistical estimations in relatively equal proportions.


In Section 2, we give a quick overview of the experiments that are typically conducted to estimate the coefficient of soilwater diffusivity. Section 3 then outlines the basic steps leading to the estimation of D(O), and Section 4 provides the details of the approach as applied to data on Manawatu sandy loam. The analysis given there also serves to illustrate (a) the exploratory use of the alternating conditional expectation (ACE) algorithm to suggest reexpressions, (b) the application of the bulge rule to aid the search for analytically tractable reexpressions of the data, (c) the use of the R 2 from the ACE transformed variables as a benchmark, and (d) the use of explicit analytic and numerical calculations to provide estimates of the derivative and integral expressions that determine D(O).
Section 5 places the transformationbased approach to the estimation of D (0) into a larger context by developing its relationship to two other problems. The first problem concerns specialized experiments for the estimation of D(O) for 0 near saturation. These experiments differ substantially from the classical horizontal infiltration experiments, and they offer an interesting direction for further statistical investigation. The second problem concerns the matching problem of reservoir engineering. It illustrates the substantial differences that can exist within the area of diffusivity estimation.


The movement of water in a horizontal column of unsaturated soil is commonly modeled by means of






the onedimensional dilfusion equation

_a2 = _ D(O) a01, at ax 1 ax
where 0 is the water content of the soil, t is time, x is the position in the horizontal column, and D(O) is the coefficient of soilwater diffusivity at the moisture level 0. Any twovariable function 0(x, t) that satisfies (2.1) can be shown (see Jost 1960, p. 31) to be a function of the single variable A = xlt"2, which is often called the Boltzman variable. After writing O(A) for the new function of one variable, one can check that O(A) satisfies the ordinary differential equation

[A] dO = d [D (0) dO]. (2.2)
2 Z Z ZA

From this equation, one can then easily obtain the 3, '4')
expression for D(O) that underlies our approach to ~(1 0 3 M/ sec' 2
its estimation: Figure 1. Scatterplot of the Volumetric Water Content, 0,
Versus the Boltzman Variable A.
D(61) dl f ),(u) du, (2.3)
2 jO determines the slope dAldO and the area under the
where 00 is the initial wcter content of the soil. The curve ~(u) for Oo s u 5 0. Naturally, there is a variety
simultaneous appearance of both derivative and in of ways in which the 0  ). plot can be used to obtain
tegral terms in this expression for D (0) provides one estimates of these functionals, but the procedure that
of the most intriguing features of its estimation. is most commonly applied is probably the one given
One point concerning (2.3) that probably deserves by Kirkham and Powers (1972, pp. 257264, 266
clarification is the absence of an additive term D(Oo) 267). A method based on splines was given by Erh
on the right side. The fact that no such term is needed (1972), and one based on piecewise parabolic fits
follows from the fact that, as 0  00, the two factors was given by DuChateau, Nofziger, Ahuja, and
of (2.3) strike an asymptotic balance. As 0  00, one Schwartzendruber (1972). Both methods provide a
has the convergence of the integral term to 0, but reasonable mechanization of the graphical determi
the derivative term goes to  in such a way that the nation of the terms of (2. 1), but neither method has
limit of the right side of (2.3) converges to D(Oo). been widely applied.
The process that has been most widely used to In addition to these socalled "graphical deter
estimate D(O) is the transientflow experiment of minations" of D(O), there are approaches that rely
Bruce and Klute (1956). In that experiment, water on estimating the parameters in an assumed para
is held at a constant head and permitted to infiltrate metric representation for D(O). One of the earliest
into a horizontal column containing airdry soil. such approaches is that of Gardner and Mayhugh
After a fixed time interval, the column is sectioned, (1958), who assumed that D(O) can be expressed in
and the water content of the individual sections is the form
determined either by weighing or by other methods.
The data of Clothier and Scotter (1982) on Mana D(O) = DdCXPOO), (2.4)
watu sandy loan plotted in Figure 1 are typical of where Dd is the soilwater diffusivity at airdry con
those obtained through horizontal infiltration exper tent, P is a parameter to be estimated, and E) is the
iments. They also give an indication of some of the dimensionless normalized water content:
inherent difficulties in estimating D(O). For instance,
many smoothing methods when applied to the data 0 = (0  00) / (0,  00), (2.5)
of Figure 1 would lead to a virtually useless estimate
of the derivative of A with respect to 0. where 0, is the saturated soil content and 00 is the
In the soilphysics literature, the determination of initial water content.
D(O) is sometimes described succinctly as follows: Other parametric models that have been explored
One plots 0 as a function of A and then "graphically" include the powerfunction form applied by Ahuja





and Schwartzendruber (1972), tion to determine the value of the definite integrals
D(O) = aOnl(O,  0)". (2.6) 1(0) = 0 A(u) du, (3.2)
This form for D(O) is suggested in part by the work
of Philip (1960), which provided a host of functional for all 00  0  o'.
forms for D(O) for which (2.1) is exactly solvable. Step 5. Let D(O) be estimated by the expression
Other parametric models were investigated by Miller [1] dl fo
and Bresler (1977), Brutsaert (1979), and Clothier D (0) =  A(u) du, (3.3)
and White (1981). Clothier, Scotter, and Green 2 dO ~
(1983) gave a method that applies for many of the where the indicated derivative and integral are those
44 exact" forms of D(O) of Philip (1960). For more determined by the results of Steps 3 and 4.

details on these methods, some comparison of their These steps are all decently explicit, except per
relative merits, and a further parametric model of haps for the first. That step does require some judg
interest, consult McBride and Horton (1985). ment, but, as the example of Section 4 illustrates,
there are useful tools to aid that judgment. Note also
3. ESTIMATION PROCESS that Step 2 takes advantage of the fact that ). is known
The description of the horizontal infiltration ex with negligible error, so we can more explicitly say
periment and formula (2.3) for D(O) enable us to that a and b are estimates in a model F(Oi) = a +
provide a topdown view of the transformation ap bG(Ai) + ei, where the error terms are independent
normals with mean 0 and variance a'.
proach to the estimation of D(O). Several features
distinguish the present approach to the estimation of 4. MANAWATU SANDY LOAM:
D(O) from previous methods, but the most salient is EXPLORATION AND TRANSFORMATION
surely the explicit role given to exploratory data anal

ysis for selecting appropriate linearizing transfor To understand the extent of the linear association
mations. The subsequent determination of the between A and 0 that can be achieved by marginal
derivative of ).(0) can then be calculated analytically, transformations, we will call on the ACE algorithm
and the integral of ).(0) can be calculated by either of Breinian and Friedman (1985). This algorithm
analytical or numerical means. finds transformations f and g such that the empirical
The details of the method we propose are possibly correlation of t',e transformed data (goi), f (0i)), 1 best explained in the context of an example such as < i:5 n, is approximately maximized. If (X, Y) is a the analysis of D(O) of Manawatu sandy loam that pair of jointly distributed random variables, one can is given in Section 4. Moreover, one almost has to define f and g as the limits of the functions f, and have an honest example to detail the role of the two g, determined by taking f 0(x) = x and go(y) = y
tools we have used to assist our transformation and applying the recursions fn 1 AX) = E(g, (Y) 1 X)
choice, the bulge rule and the ACE algorithm. With and g,,,(Y) = E(f,,(X) 1 Y). The f and g deter
that said, it seems useful to have a topdown view mined by this process can be shown to maximize
of the proposed method. The five basic steps are the corr(f (X), g(Y)). Naturally, if the joint distribution
following: of X and Y is not known, one cannot find f and g
Step 1. Find transformations F(O) and Go.) such precisely by this method, but one can derive an em
that the transformed data values (F(Oi), G(Ai)), 1  pirically based algorithm as a natural modification of
i  n, exhibit a strong linear association. This step the theoretical algorithm. All of our ACE compu
will be achieved by means of the ACE algorithm and tations were performed usingthe implementation of
the empirical ACE algorithm of L. Breiman that was
the bulge rule for selecting power transformations. incorporated by Schilling (1985).

Step 2. Use simple linear regression (or a more Even for the best choices of f and g, the linear
sophisticated technique like iteratively reweighted association between ). and 0 is imperfect. Still, the
least squares) to determine values a and b, which ACEtransformed variables plotted in Figure 2 ex
provide an approximate functional relationship in 0 hibit substantially greater linear association than the
and ). of the form plot of the untransformed variables given in Figure
F(O) = a + bGo.). (3.1) 1. Moreover, when we measure the linear association
of the transformed variables in terms of R', we find
Step 3. Use a chain rule calculation to extract from a respectable value of R' = .93. This value provides
(3.1) an expression for dAldO in terms of 0. us with a benchmark; in fact, one of the principal
Step 4. Use either analytic or numerical integra benefits of the ACE algorithm is that it provides a


............ ~ ~


2 1 0

.. T

Figure2. Scatterplot of the A CE Transformed Oi Versus the ACETransformed Ai. This plot is used to assess linearity of the ACE transformation.

theoretical standard against which more analytically appealing transformations can be judged.
To aid in the search for such surrogates for f and g, the ACEtransformed variables are plotted against the untransformed variables to see if simpler functional forms might suffice. Figures 3 and 4 show the


1 1 1

0 1 2



4 5


Figure 3. Scatterplot of the ACETransformed Ai Versus Ai That Is Used to Suggest a Power Transformation Approximating g (A).


i i

i i

zz:~ 7 


1 1 1 1

0.10 0.20 0,30


Figure 4. Scatterplot of the ACETransformed Oi Versus Oi That Is Used to Suggest a Power Transformation Approximating f (0). Notice that the bulge rule may not be directly applicable here.

plots of (Ai, g(Ai)), 1 :5 i  n, and (0j, f (0)), 1  i :5 n, respectively.
The hunt for analytically tractable replacements for f and g is further guided by the socalled bulge rule of Mosteller and Tukey (1977). Loosely speaking, the bulge rule suggests finding an outward normal to a smoothed plot of the data and using the signs of the normal components to guide one's choice of transformation. For example, Figure 3 exhibits a bulge where both the x and y components of the outward normal are positive. The bulge rule then suggests that the variable plotted on the horizontal axis should be transformed by moving up the scale of powers. In fact, the successive examination of plots of (A7, g(Ai)), 1 :5 i :5 n, for larger values of a continues to suggest moving up the scale of powers, and we are thus led to consider the exponential transformation. An alternative approach to this exploratory search for an appropriate transformation would be to use the method of Box and Cox (1964). The results obtained by the BoxCox method are comparable to those obtained via the bulge rule.
When we begin a similar examination of the plot of (0j, f (0)), 1  i :5 n, given in Figure 4, the bulge rule for reexpression fails to offer unambiguous advice. There may be a modest indication that we might wish to send 0 down the scale of powers, but the indication is not supported when tried. Fortunately, we have recourse to a second approach that does suggest an appropriate transformation, and we can



0 100 200 300 400


Figure 5. Scatterplot of Oi Versus e,.i. The bulge rule suggests going up the ladder for either 0 or A.

consider the plot of Oi versus e~,, which is shown in Figure 5. After all, since we have settled on eA as the surrogate for f (~), the principal remaining task is to determine a surrogate F for f such that the scatterplot (1'(0), e~i) is approximately linearized.
The bulge rule applied to Figure 5 initially suggests that we c,,nsider a transformation F that moves 0 up


the ladder of powers, and successive applications of the bulge rule eventually lead us to the choice of F(O) = 03. Strikingly, this transformation achieves an R` of .93 that meets the level of the optimal R' = .93 achieved by the ACE transformations. Moreover, when we consider the plot of 0,~ versus eli given in Figure 6, the visual impact of the linear association exhibited by this figure seems to compare well with that exhibited by the ACEtransformed variables of Figure 2. On the basis of both the quantitative evidence provided by comparing R 21 s and the subjective evidence provided by comparison of the scatterplots of Figures 2 and 6, it seems appropriate to settle on the transformation choices of Figure 6.


For the Manawatusandyloam data, our exploratory analysis has led us to an approximate relationship of the form

F(O) = a + bG(A), (5.1)

where F(O) = 03 and G(A) = el. The coefficients in (5. 1) can now be estimated by ordinary least squares, from which we obtain a = 4.48 x 10~1 and b =  1.20 X 104 with nominal standard errors of 5.30 X 104 and 3.30 x 106, respectively. In Figure 7 we show a plot of the predicted values O~ versus the residuals that are obtained from fitting the model (5. 1) by ordinary least squares. The residuals appear approximately homescedastic, and we have no rea






Figure 6. Scatterplot of 0,1 Versus ei. Approximate linearity is achieved with this transformation, which should be compared with the ACE transformations of Figure 2.



12 q

1 1 1 1 1
0 100 200 300 400 0.0 0.01 0.02 0,03 0.04

Predicted 0

Figure 7. Scatterplot of Residuals Versus Predicted Values for the Model 07 = a + beli. Residuals appear to be approximately homoscedastic.




son to be discontent with the estimates obtained by ordinary least squares. If the scatterplot of Figure 7 had exhibited a greater heteroseedasticity, we would have probably elected to apply iteratively reweighted least squares or a similarly directed technique.
As a final check on the reasonability of the fitted model, one should consider the fit in terms of the untransformed variables as exhibited in Figure 8. This plot has no flagrant defects; indeed, it suggests that the procedure has been reasonable. Still, one feature of Figure 8 deserves comment. The principles of good experimental design suggest that one should take more observations in those regions of A where 0 changes more rapidly. Practitioners are aware of this fact and often take thinner soil sections near the wet front, although no systematic analysis of the experimental design has yet been provided.
By differentiating (5.1), we find the general relationship

dA = F'(0) = b1F'(0) ~ (5.2)
dO W(A) W(G1QF(O)  a)lb))

and, for the choices that were made by means of the exploratory analysis of the Manawatusandyloam data, one finds a particularly simple net result:

dA = 02(01  j)1.
j 3 (5.3)

To obtain D(O), it remains only to determine the integral of A(O) = G1QF(O)  a)lb). For the Manawatusandyloam data, we find A(O) = logQ01 




g ~C r p 5"


0.10 0.15 0,20
0 1 2 3 4 5 6 0

Figure 8. Scatterplot of Oi Versus A, Predicted values from the equation & = (a + be4)113 are shown by the line.




2 1 ,~? 6




1 1 1
0.10 0.15 0.20


 7   1 1 1
0.25 0.30 0.35

Figure 9. Plot of SoilWater Diffusivity D(O) Versus Moisture Content 0.

a)lb), and the integral of A(O) can be determined analytically.
We are now in the position to provide the estimates for D(O). These are given first in Figure 9 on a raw scale and again in the more appropriate log scale in Figure 10. The 99% confidence intervals given in Figure 10 probably deserve some comment. By the

1 1 r_
0.25 0,30 0,35

Figure 10. Plot of Loglo of SoilWater Diffusivity Versus 0 and the Associated 99% Upper and Lower Confidence Intervals.


usual least squares theory, one can obtain the joint recognized by Bruce and Klute (1956) and were distribution of the estimates j and 6 of the param further emphasized by MorelSeytoux and Khanji eters of Equation (5. 1). We then calculated 10' pseu (1975). More recently, experiments to address the dorandom observations (ai, bi) with the same joint estimation of D(O) near saturation were performed distribution as (d, 6). For each 0 in a set of 100 evenly by Clothier and Wooding (1983) and applied by spaced values, we then calculated Di(O) based on Clothier et al. (1983). These new experiments are (ai, bi). Finally, we set D*(0) and D*(0) equal to substantially more sophisticated than the classical the upper and lower 1% points of {Di(O): 1 :5 i .5 horizontal infiltration experiments. They rely on the 10% and used these values to provide our estimates use of continuous dripping of water onto a column of the upper and lower values of the confidence of soil and observation of the associated periodic bounds in Figure 10. The extreme narrowness of fluctuations in pressure potential at differing depths these confidence bounds is fundament lly a reflection in a vertical soil column. These procedures appear
of the small standard errors of d and b. One essential promising, and reexamination of the associated es
observation concerning these confidence bounds is timations from a statistical viewpoint would be val
that their width increases greatly as 0 approaches the uable.
saturation level 0,. This feature is endemic to the When one looks beyond the immediate problem
problem of estimating D(O), and one would have of estimating the coefficient of soilwater diffusivity,
been alarmed if the confidence bounds on D(O) did one finds a world of estimation problems associated
not degrade as 0 approached saturation. with diffusion models. Some of these problems are
The desire to know how well we have done is close cognates to the problem studied here, but other
natural but frustrating. In fact, one faces today al problems with a strong structural resemblance are
most the same situation as when Bruce and Klute faux amis. A striking example of this phenomenon
(1956) observed, "There is unfortunately no stan is provided by the historymatching problem of res
dard against which diffusivity values can be checked" ervoir engineering that was discussed from a statis
(p. 462). The estimates of Figures 9 and 10 pass the tical point of view by O'Sullivan (1986).
tests of physical reasonability and exhibit reasonable Again, the basic problem concerns the estimation
consistency with the results obtained by other meth of the diffusion parameter in a model of flow through
ods, such as those of McBride and Horton (1985). a porous medium. Specifically, one wants to estimate
We are, therefore, at a point that has confronted all the coefficient a(x) in the diffusion equation with
others who have attempted to estimate D(O). We forcing:
have an estimate that appeals to our theoretical sen a U (X, t) a au
sitivities, but, on the basis of current experimental a (x)  (X, t) q (x, t),
evidence, we cannot argue that the present estimate at ax ax
is indisputably better than those obtained earlier. x EE 1, 0 :s t :s T. (6. 1)
This fact is part of the interactive process between
theory and experiment, and from the beginning of Here x is a twodimensional spatial variable, and one
this investigation no outcome could have been ex has data on scattered well pressures u(xi, t) and
pected to emerge with an ironclad claim on the truth. forcing pressures q(xi, b), 1 :5 i :5 m, 1 25 j :5 1.
6. CONCLUDING OBSERVATIONS Three essential facts separate the reservoirmatch
ing problem from the problem of estimating soil
The problem of estimating the coefficient D(O) of water diffusivity. First, the spatial variable is necessoilwater diffusivity is one that deserves the atten sarily two dimensional, so there is no possibility of tion of statisticians. The approach explored here is a reduction to an ordinary differential equation like natural from the point of view of contemporary data (2.2). Second, time plays an essential role in (6.1), analysis; yet in comparison to commonly applied whereas in the analysis of (2.1) time had only the
analyses, the present approach offers some substan modest role of scaling the Boltzman variable A =
tial benefits. First, the estimate for D(O) is guaran x1P`. Finally, in (6. 1) one has to deal with a complex
teed to be smooth and monotone, just as one would forcing term. The presence of a modest forcing term
expect from fundamental physical principles. Fur does not automatically remove one from the domain
ther, the present method permits one to be genuinely of the present analysesfor example, the methods
guided by the data in contrast to methods that assume applied here can be modified for use on vertical in
parametric forms for D(O) based on past experience filtration experiments in which one has a forcing term
or analytic convenience. because of gravity. But forcing terms with the com
Still, it should be understood that the estimation plexity of those present in the reservoirmatching
of D(O) is not always easy, especially for values of problem are of a different order, and their presence
0 near saturation. The difficulties in that range were takes one into a different realm.



The two problems just sketched should serve to tration From a Hemispherical Cavity," Soil Science Society of
suggest that the subject of diffusivity estimation is America Journal, 46, 696700.
quite broad. In fact, our original problem of esti Clothier, B. E., Scotter, D. R., and Green, A. E. (1983),Dif
fusivity and One Dimensional Absorption Experiments," Soil
mating soilwater diffusivity lies at the beginning of Science Society of America Journal, 47, 641644.
a whole scale of challenging problems that appear Clothier, B. E., and White, 1. (1981), Measurement of Sorptivity
ripe for statistical exploration. We hope that the steps and SoilWater Diffusivity in the Field," Soil Science Society of
taken here provide some introduction and encour America Journal, 45, 241245.
agement to statisticians who would engage the many Clothier, B. E., and Wooding, R. A. (1983), "The SoilWater
Diffusivity Near Saturation," Soil Science Society of America
difficult diffusivity estimation problems that remain. Journal, 47, 636640.
DuChateau, R C., Noffiger, D. L_ Ahuja, L. R., and Schwartz
ACKNOWLEDGMENTS endruber, D. (1972), "Experimental Curves and Rates of
Change From Piecewise Parabolic Fits," Agronomy Journal,
We thank John MeBride for his help in introducing 64, 538542.
us to this problem and for providing the data sets Erh, K. T. (1972), "Application of the Spline Function to Soil
Science," Soil Science, 114, 333338.
that initiated our analysis. We also thank B. E. Cloth Gardner, W. R., and Mayhugh, M. S. (1958), "Solutions and Tests
ier for his kind comments on an earlier draft. This of the Diffusion Equation for the Movement of Water in Soils,"
research was supported in part by National Science Soil Science Society of America Proceedings, 22, 197201.
Foundation Grant DMS8414069. Jost, W. (1960), Diffusion in Solids, Liquids, and Gases (3rd ed.),
New York: Academic Press.
[Received March 1988. Revised May 1988.1 Kirkham, D., and Powers, W. L. (1972), Advanced Soil Physics,
New York: WileyInterscience.
McBride, J. E, and Horton, R. (1985), "An Empirical Function
REFERENCES to Describe Measured Water Distributions From Horizontal
Infiltration Experiments," Water Resources Research, 21, 1539
Ahuja, L. R_ and Schwartzendruber, D. (1972), "An Improved 1544.
Form of SoilWater Diffusivity Function," Soil Science Society Miller, R. D., and Bresler, E. (1977), "A Quick Method for
of America Proceedings, 36, 914. Estimating SoilWater Diffusivity Functions," Soil Science So
Box, G. E. R, and Cox, D. R. (1964), "A., Analysis of Trans ciety of America Journal, 41, 10201022.
formations" (with discussion), Journal of the Royal Statistical MorelSeytoux, H. J_ and Khanji, J. (1975), "Pred:,~tion of 'in
Society, Ser. B, 26, 211252. hibition in a Horizontal Column," Soil Science Society ofAmer
Breiman, L_ and Friedman, J. (1985), "Estimating Optimal ica Proceedings, 39, 613617.
Transformations for Multiple Regression and Correlation," Mosteller, E, and Tukey, J. W. (1977), Data Analysis and Regres
Journal of the American Statistical Association, 80, 580597. sion, Reading, MA: AddisonWesley.
Bruce, R. R_ and Klute, A. (1956), "The Measurement of Soil O'Sullivan, F. (1986), "A Statistical Perspective on 111Posed In
Moisture Dilfusivity," Soil Science Society of America Proceed verse Problems," Statistical Science, 1, 502527.
ings, 20, 458462. Philip, J. R. (1960), "General Method of Exact Solution of the
Brutsaert, W. (1979), "Universal Constants for Scaling the Ex Concentration Dependent Diffusion Equation," Australian
ponential Soil Water Diffusivity?" Water Resources Research, Journal of Physics, 13, 112.
15, 481483. Schilling, J. M. (1985), The Statistics Store, Murray Hill, NJ:
Clothier, B. E_ and Scotter, D. R. (1982), "ConstantFlux Infil AT&T Bell Laboratories.