## Report DOE/ER/82329-2, The Brooks-Corey Approximation

prev   index   next

Notes: This section, as an excerpt, refers to figures not presented on this web site. You must order the full publication to see them. It will be helpful to go to Gathered Equations and print out the equations before reading this section. This page shows greek characters correctly with Netscape Navigator 3.0 on a Win 95 computer with the SYMBOL TrueType font set. If it does not do so on your computer, then psi may show as y, eta as h, and delta as D.

Approximate Darcian Means for Brooks-Corey Conductivity Relations

Darcian means get much more interesting for a Brooks-Corey conductivity relation [34]. With an exponential relation, there's not much to plot because the relevant weighting functions, wxo [42] and wx1 [29], do not change over the entire (kr1,kr2)-space. To use the weighting of [38] for a Brooks-Corey relation, we can develop an analytic expression of Kh [45] by substituting [34] into [8]. We will call this special function, km, the "manta function", due to its shape in log3 space. It has a number of special cases, shown in Table 5, which include most standard means. It provides a smooth transition from the harmonic mean, m=-2 to the arithmetic mean, m=1, and illustrates that most standard means are only special cases in an infinite continuum of means.

[45]

Table 5: Some Standard Means as Special Cases of the Manta Function

 m km standard mean 1 arithmetic ½ … 0 log -½ geometric -1 … -2 harmonic

The analytic derivation of an equivalent expression for kxo [41] in this case is not so easy. The distribution, kr(x), for a Brooks-Corey relation in this context seems at best to be an implicit solution, if any in fact exists. This investigator made plots of kx(kr1,kr2) for Brooks-Corey relations [34] with h varying from 1 to 46 and kr varying from 10-10 to 1, for Dx very small, from solutions of [11] for which the relative extreme flow difference [17] was good to about 1 part in 106. It turns out that km(m=1-(1/h)) fits kx in most of these cases to a relative error [46] of about 1 part in 104. This formula for relative error has been very useful in the huge dynamic range of this study.

[46]

Coincidentally, one can also get the same expression for km by assuming that kr(x) = (a+bx)1/m and substituting into [7]. But it is hard to find the reasoning that would make this an analytic derivation. The reader who can demonstrate it should publish. This does, however, show that all such means based on this assumption are valid only for Dx à 0.

Let us see what we can determine by plotting Kh, Kv and Kx for relation [34] for the Brooks-Corey parameters for separate fracture flow (yd = 0.082 m, h = 8.075) and matrix flow (yd = 33m, h = 1.9825) in a Topopah Spring welded tuff. Figures 19-23 were developed in much the same way as Figures 3 and 4 from numerical solutions to the elliptic boundary condition problem [11]. Figures 19a & b show Kh for the fracture and matrix flow, respectively. Figures 20a-e show the total Darcian mean, Kv, for the fracture flow for decade values of Dx at 0.00082 to 8.2m. Figures 21a-e show Kv for the matrix flow for Dx at decades from 3.3 to 3300m. Figures 22a-e and 23a-e show the Darcian gravity flow mean, Kx, for the same conditions of fracture and matrix flow, respectively. The reader will notice that Figures 20d, 20e and 22d have outliers. The effect on the surface fitting-software used to plot these data was so bad in 20d that the mesh surface was left out.

 Figure 19a: Kh for fracture flow Figure 19b: Kh for matrix flow

 Figure 20a: Fracture flow Kv for Dx = 0.00082m Figure 21a: Matrix flow Kv for Dx =0.33m Figure 20b: Fracture flow Kv for Dx = 0.0082m Figure 21b: Matrix flow Kv for Dx = 3.3m Figure 20c: Fracture flow Kv for Dx = 0.082m Figure 21c: Matrix flow Kv for Dx = 33m Figure 20d: Fracture flow Kv for Dx = 0.82m Figure 21d: Matrix flow Kv for Dx = 330m Figure 20e: Fracture flow Kv for Dx = 8.2m Figure 21e: Matrix flow Kv for Dx = 3300m

 Figure 22a: Fracture flow Kx for Dx = 0.00082m Figure 23a: Matrix flow Kx for Dx = 0.33m Figure 22b: Fracture flow Kx for Dx = 0.0082m Figure 23b: Matrix flow Kx for Dx = 3.3m Figure 22c: Fracture flow Kx for Dx = 0.082m Figure 23c: Matrix flow Kx for Dx = 33m Figure 22d: Fracture flow Kx for Dx = 0.82m Figure 23d: Matrix flow Kx for Dx = 330m Figure 22e: Fracture flow Kx for Dx = 8.2m Figure 23e: Matrix flow Kx for Dx = 3300m

Figures 19a, 20a, 22a and 23a demonstrate the reason for the term "manta function". The plot looks a bit like a manta ray with its wings up, and the spine laying along the kr1=kr2 diagonal. In Figure 19b, with h = 1.9825, the value of m = -1/1.9825 is very close to that for a geometric mean (Table 5) and thus the wings are characteristically flat. In plots 20, 22 and 23, m is significantly greater than -1/2, about -0.12384, 0.87616 and 0.49559, respectively, so that the wings are always up. For a harmonic mean, the wings would be down. Because for kxo, m = 1-1/h > -1/h, the curl at the spine is sharper and the wings rise higher for Kx than for Kv.

Note that Figures 20a and 21a are practically indistinguishable from Figures 19a and b, as Dx is small enough to make Kv look much like Kh. Likewise, Figures 22a and 23a are close to what kxo would look like for the fracture and matrix flows. Note how, with the exception of Figures 21a-e, as Dx increases over orders of magnitude, the spine of the manta tends to roll to the right, into the region of dry-over-wet. At the same time, the left wing rises into the kr2 position, demonstrating the transition from the prototype forms, Kh and kxo, for Dx à 0 to the up-gravity, or kr2 mean. Again, Kv seems to make this transition well ahead of Kx.

It may seem silly to consider values of Dx ranging from 820 microns to 3.3 kilometers, but this illustrates how Darcian means change for different materials at different model scales. For example, consider the Darcian gravity flow means, Kx, in Figures 23b, 22c and 22d. At Dx = 0.3 m, it would seem that kx for the matrix flow (Fig. 23b) is relatively close to the manta function limit, km(m=0.496), close to (kr1+kr2+(kr1*kr2)1/2)/3 in Table 5. In Figures 22c and d, Kx for fracture flow on the left hand side, wet-over-dry conditions, looks to be close to the kr2 limit. On the other side of the kr1-kr2 diagonal, for dry-over-wet conditions, the curl of the spine precludes the kr2 mean for Kx.

While these figures are useful in understanding how Darcian means evolve with vertical space step size, they are not as useful in constructing adequate approximations. Consider Figure 22e. It looks almost as though the spine of the manta is parallel to the kr1-kr2 diagonal. It's not. It's about 3 to 4 degrees off. This makes it practically impossible to find an analytic transformation that will take Kx from kxo to kr2 by any kind of rotation or scaling in log3-space. The same is true for Kv. Besides that, Figure 21d demonstrates a saddle shape that kills any prospect of that kind of transformation

But, as it happens, if we plot Kx(kr1,kr2) and Kv(kr1,kr2) not for constant-Dx surfaces, but for constant-u surfaces (see [48]), everything lines up parallel to the kr1=kr2 diagonal. This allows us to do something very useful. For a single set of parameters, (yd, h, u), we can describe Kx or Kv each as a single curve in a new 2-D log space that looks down the kr1=kr2 diagonal. We can collapse every 3-D surface into a 2-D curve.

Again, the reader who can demonstrate this property through analytic derivation is welcome to publish on it. It's not an easy problem. We are talking about parameters and relationships developed from literally thousands of solutions to very nonlinear elliptic boundary condition problems that have no readily apparent analytic solution themselves. And those who wish to comment on not proving all the cases have only to solve such problems and demonstrate the exceptions.

So, how can we look down a tilted diagonal in log3-space and transform a surface to a curve? Consider the transformation in [47] from the (kr1,kr2)-plane in log3-space to the (tl,tw)-plane. The long direction, tl, is parallel to the kr1=kr2 diagonal. The wide direction, tw, measures relative distance in log3-space away from the diagonal. The directions of tl and tw are normal in log3-space and cause the desired collapse from 3-D constant-u surfaces to 2-D curves.

[47]

But what does constant-u mean? Equation [36] has only one value of y, and we have both y1 and y2 at x1 and x2. It turns out that we can convert to kr through equation [35] and substitute tl for kr as shown in [48] and [49]. In effect, we use the geometric mean of kr1 and kr2. This allows the collapse of 3-D surfaces to 2-D curves.

[48]

[49]

Figures 24a-d show log(Kx/tl) plotted against log(kr2/kr1) (or log(tw2)), with u varying in 51 log steps from 0.025 through 3.593512 to 516.53. Figures 25a-d show Kv in the same conditions. Each of these curves is a tilted plane in log3-space. The kr1-kr2 diagonal is the end-on point in the center of the plots. Note that in each of these plots, Kx and Kv transform from the prototypical kxo or Kh to kr2, the diagonal line from lower left to upper right. The constant-u curves spread out in the middle and bunch toward the limits because of the nonlinear behavior of [42] and [43]. Although the analytic connection cannot as yet be demonstrated, there is still an effect.

A few notes: The Kx plots curve more sharply upward than the Kv plots, because for the same h, mx=mv+1. The curves plotted more smoothly in the original spreadsheet; the Window 95 clipboard cut and paste is not perfect and introduces some jaggies. See how all the curves pass through the kr1=kr2 diagonal, the point in the center of the plots, a necessary condition for a valid interblock conductivity mean. Note that except through the effect of the scaling parameter, u, on Dx, the "displacement pressure", yd, does not affect these curves.

 Figure 24a: Kx/tl for h = 2/3, m = -1/2 Figure 25a: Kv/tl for h = 2/3, m = -3/2 Figure 24b: Kx/tl for h = 1.9825, m = 0.4956 Figure 25b: Kv/tl for h = 1.9825, m = -0.5044 Figure 24c: Kx/tl for h = 8.075, m = 0.8762 Figure 25c: Kv/tl for h = 8.075, m = -0.1238 Figure 24d: Kx/tl for h = 20, m = 0.95 Figure 25d: Kv/tl for h = 20, m = -0.05

This approach also gives us a more comprehensive way to visualize the error of standard means and approximations to Darcian means. For example, consider the Darcian means in Figures 24c, Kx, and 25c, Kv, for a Brooks-Corey relation with h = 8.075 for fracture flow in a Topopah Spring welded tuff. Figures 26a-d show the ratio of means ka, kg, kxo and kr2, respectively to Kx. Figures 27a-d show the ratios of means ka, kg, Kh and kr2 to Kv. Kh corresponds to the integral mean, ki, in the beginning of this work, ka is the arithmetic mean, kg is the geometric mean, and the upstream mean is neglected. The plots show curves of constant-u for u varying from 0.025 to 516.53 in 21 log steps.

 Figure 26a: Ka error for h = 8.075, Kx Figure 27a: Ka error for h = 8.075, Kv Figure 26b: Kg error for h = 8.075, Kx Figure 27b: Kg error for h = 8.075, Kv Figure 26c: Kxo error for h = 8.075, Kx Figure 27c: Kh error for h = 8.075, Kv Figure 26d: Kr2 error for h = 8.075, Kx Figure 27d: Kr2 error for h = 8.075, Kv

For the given range of kr2/kr1, the mean ka, kg, kxo, Kh, or kr2, matches Kx or Kv only on the horizontal log = 0 line. Everything off that line is in error. Note that this occurs mathematically only in Figures 26c and 27c, and comes quite close in Figures 26a and 27d, and only occurs at one end of the range of u. The (0,0) center in these log-log plots is the kr1=kr2 diagonal in the 3-D plots exactly. Note how fast many of the means diverge from ideal as one moves away from the kr1=kr2 diagonal. This does not bode well for models using these means.

Now let us consider whether or not weighted approximations [27] and [28] using the weighting definition [29], which were exact for the case of an exponential conductivity relation, work as well for the case of a Brooks-Corey relation. Let us also consider the approximations [38] and [39], using the weighting definition [42]. In all cases, we are using the definition [48] to develop u from (kr1,kr2).

Let kxw1 be the value of Kx estimated from [27], kvw1 the value of Kv estimated from [28] for wv1 = wx1, kxwo the value of Kx estimated from [39] and kvwo the value of Kv estimated from [38] for wvo = wxo. Figures 28a-d show the ratio of kxw1/Kx versus kr2/kr1 in log-log format, Figures 29a-d the ratio of kvw1/Kv, Figures 30a-d the ratio of kxwo/Kx, and Figures 31a-d the ratio of kvwo/Kv, for h = 2/3, 1.9825, 8.075 and 20, respectively for a-d. The values for Kx and Kv are taken from the data for Figures 24 and 25. Because there is no real difference in kxw1 and kvw1, the only difference between Figures 28 and 29 are Kx and Kv used in the ratios. The real difference in kxwo and kvwo is the prototypical Kxo and Kh in [39] and [38].

 Figure 28a: kxw1/Kx for h = 2/3 Figure 29a: kvw1/Kv for h = 2/3 Figure 28b: kxw1/Kx for h = 1.9825 Figure 29b: kvw1/Kv for h = 1.9825 Figure 28c: kxw1/Kx for h = 8.075 Figure 29c: kvw1/Kv for h = 8.075 Figure 28d: kxw1/Kx for h = 20 Figure 29d: kvw1/Kv for h = 20

 Figure 30a: kxwo/Kx for h = 2/3 Figure 31a: kvwo/Kv for h = 2/3 Figure 30b: kxwo/Kx for h = 1.9825 Figure 31b: kvwo/Kv for h = 1.9825 Figure 30c: kxwo/Kx for h = 8.075 Figure 31c: kvwo/Kv for h = 8.075 Figure 30d: kxwo/Kx for h = 20 Figure 31d: kvwo/Kv for h = 20

Note how the curves for kxw1, kvw1, kxwo and kvwo, for h = 8.075, compared to the true Darcian means, stay much closer to the ideal ratio of 1:1 farther from the kr1=kr2 condition than do the standard means in Figures 26 and 27. Note that since we have not attempted, in this case, to account for the variation of wv1 and wvo away from the kr1=kr2 diagonal, as indicated in [30], the estimates for Kv tend to be worse than the estimates for Kx. Note how for smaller values of h, kxwo tends to be better than kxw1, and for all the values of h shown, kvwo tends to be better than kvw1.

So, this provides us with a reasonable set of first approximations to Darcian means for the Brooks-Corey relation. They should allow the model to move farther in state space from the kr1=kr2 diagonal with less error in calculated flow. Since the Brooks-Corey relation [34] is determined by two parameters, then like the exponential relation [22], it, too, can be used in piecewise manner in a model to estimate the Darcian means between adjacent grid points for more complicated relations. If we set up the piecewise approximation in Figure 1, for [34] as we did for [22], we get h [50] and yd [51].

[50]

[51]