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.
Darcian means get much more interesting for a BrooksCorey 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 BrooksCorey 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 log^{3} 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]





















The analytic derivation of an equivalent expression for kxo [41] in this case is not so easy. The distribution, kr(x), for a BrooksCorey 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 BrooksCorey 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 10^{6}. 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 10^{4}. 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 BrooksCorey 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 1923 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 20ae show the total Darcian mean, Kv, for the fracture flow for decade values of Dx at 0.00082 to 8.2m. Figures 21ae show Kv for the matrix flow for Dx at decades from 3.3 to 3300m. Figures 22ae and 23ae 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 fittingsoftware used to plot these data was so bad in 20d that the mesh surface was left out.


























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 = 11/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 21ae, as Dx increases over orders of magnitude, the spine of the manta tends to roll to the right, into the region of dryoverwet. 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 upgravity, 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, wetoverdry conditions, looks to be close to the kr2 limit. On the other side of the kr1kr2 diagonal, for dryoverwet 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 kr1kr2 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 log^{3}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 constantDx surfaces, but for constantu 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 2D log space that looks down the kr1=kr2 diagonal. We can collapse every 3D surface into a 2D 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 log^{3}space and transform a surface to a curve? Consider the transformation in [47] from the (kr1,kr2)plane in log^{3}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 log^{3}space away from the diagonal. The directions of tl and tw are normal in log^{3}space and cause the desired collapse from 3D constantu surfaces to 2D curves.
[47]
But what does constantu 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 3D surfaces to 2D curves.
[48]
[49]
Figures 24ad show log(Kx/tl) plotted against log(kr2/kr1) (or log(tw^{2})), with u varying in 51 log steps from 0.025 through 3.593512 to 516.53. Figures 25ad show Kv in the same conditions. Each of these curves is a tilted plane in log^{3}space. The kr1kr2 diagonal is the endon 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 constantu 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, m_{x}=m_{v}+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.








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 BrooksCorey relation with h = 8.075 for fracture flow in a Topopah Spring welded tuff. Figures 26ad show the ratio of means ka, kg, kxo and kr2, respectively to Kx. Figures 27ad 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 constantu for u varying from 0.025 to 516.53 in 21 log steps.








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 loglog plots is the kr1=kr2 diagonal in the 3D 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 BrooksCorey 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 28ad show the ratio of kxw1/Kx versus kr2/kr1 in loglog format, Figures 29ad the ratio of kvw1/Kv, Figures 30ad the ratio of kxwo/Kx, and Figures 31ad the ratio of kvwo/Kv, for h = 2/3, 1.9825, 8.075 and 20, respectively for ad. 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].
















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 BrooksCorey 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 BrooksCorey 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]