A simple method for solving the Bussian equation for electrical conduction in rocks

Abstract. One of the most general and effective models for calculating the complex electrical conductivity and relative dielectric permittivity of rocks saturated with pore fluids is that of Bussian. Unlike most models, it is non-linear and cannot be solved algebraically. Consequently, researchers use reiterating numerical routines to obtain a solution of the equation, and then only for the real part of the solution. Here we present a different approach to the solution that uses conformal mapping in the complex plane, and implements it within MapleTM. The method is simple and elegant in that it requires, for example, only 3 lines of code in MapleTM 11 and little programming experience. The approach has been shown to be as precise as using the classical reiterating bisection method for real data implemented in C++ on an ordinary desktop computer to within a probability over 1 in 109. However, the conformal mapping approach is 52 times as fast. We show once more that the Bussian equation breaks down for low fluid conductivities, but recommend it (with the modified Archie's law) for use with rocks saturated with high salinity fluids when the matrix is conductive.


Introduction
The measurement and understanding of the electrical conductivity of porous media has applications in many areas of science and technology.Perhaps the most important is its use in the oil and gas industry for calculating the reserves of hydrocarbons in reservoirs from electrical well logging measurements.The original relationship for interpreting resistivity logs is Archie's law (Archie, 1942), which was arrived at empirically.Its use is restricted to rocks where water saturating the pores is the only conductor (i.e., clean sandstones Correspondence to: P. W. J. Glover (paglover@ggl.ulaval.ca)and carbonates without conductive accessory minerals at low temperatures).Archie's law has been modified to extend its range of applications many times: (i) for a matrix of water saturated conductive particles (Wyllie and Southwick, 1954;de Witte, 1957).Later, Waxman and Smits (1968) proposed the first model for shaly sand formations.In 1977, Clavier et al. (1977) suggested a model that assumes that Stern layer cations contribute to the conductivity of the clay water and bulk water separately, and is consequently called the dualwater model.More recently Archie's law has been extended to work for two conducting phases (Glover et al., 2000a), and is particularly useful for applications at high temperatures where the matrix has a finite conductivity (Glover et al., 2000b).Even more recently a generalised Archie's law for n phases has been published (Glover, 2010).
Until 1980 Archie's law had an empirical pedigree.Then it was shown by Sen (1980Sen ( , 1981) ) and Sen et al. (1981) that Archie's law follows from the work of Bruggeman (1935) and Hanai (1960aHanai ( , b, 1961)), which is a well-defined theoretical model that is based upon classical physics and geometry and assumes that non-conducting rock particles are dispersed in a continuous phase of saline water.The Bruggeman-Hanai-Sen (BHS) equation can be expressed as where the effective complex conductivity of the mixture σ * eff is expressed relative to the complex electrical conductivity of dispersed particles σ * p within a continuous medium with a complex electrical conductivity σ * f .Here φ p is the volume fraction of particles and d is the so-called depolarisation factor.The conductivities are complex and follow the relationship σ * = σ + iωε o κ, where ω is the angular frequency, ε o is the electric constant (ε o = 8.854 × 10 −12 F/m), κ is the relative dielectric constant, and i = √ −1.Hence the equation implicitly includes AC current transport.
Published by Copernicus Publications on behalf of the European Geosciences Union.P. W. J. Glover et al.: A simple method for solving the Bussian equation for electrical conduction in rocks Bussian (1983) reinterpreted the BHS approach to include a conducting lattice-like matrix.His model relates the electrical properties of any heterogeneous two-component mixture to the properties of the individual components at any frequency, and replaces the depolarisation factor d with a parameter m = 1/(1d), that can be shown to be the same as the classical Archie cementation exponent.In the Bussian equation the complex effective conductivity σ * eff , the complex dielectric permittivity ε * eff , and the complex relative permittivity κ * eff of a two phase medium follow the equations where m is the cementation exponent, which describes the effect that the arrangement of the pore space has on the electrical parameters, σ * m , ε * m and κ * m are the complex conductivity, complex dielectric permittivity and complex relative permittivity of the matrix, σ * f , ε * f and κ * f are the complex conductivity, complex dielectric permittivity and complex relative permittivity of the fluid that saturates the pores, and φ is the porosity.
This equation is (i) general, treating the electrical transport properties as complex parameters, (ii) valid for all frequencies, and (iii) reduces to classical laws for special cases.It should be the equation of choice when modelling the electrical properties of porous media with saline fluids.However, the equation is not valid at low fluid conductivities.This limitation arises from its origins in effective medium theory (Bussian, 1983).Unfortunately, the equation is non-linear.The definition of a non-linear system is one in which the variable(s) to be solved for cannot be written as a linear combination of independent components, i.e., it is a system which does not satisfy the superposition principle.Hence, the Bussian equation cannot be solved algebraically.The complications and limitations involved in solving the equation using reiterative numerical methods means that it is often overlooked.
Another method that is based on the mixing of conductivities of two components, and that is derivable analytically, is that of Korvin (1982) and Tenchov (1998).We later compare the Bussian equation with the modified Archie's law and the Korvin and Tenchov approach.

Conventional approaches to solving Bussian's equation
There are many methods available for the solution of nonlinear equations numerically and it is beyond the scope of this paper to review them all.An excellent and accessible review of all of the methods discussed below is available together with code in Numerical Recipes in C++: The Art of Scientific Computing, 3rd edition (Press et al., 2007), in which the bisection method is described as an extremely robust method that cannot fail for smoothly varying well defined functions.They also recommend the Brent method and Ridder's method especially if the function cannot easily be differentiated, as in our case.Differentiable functions can make use of the Newton-Raphson method with Press et al.'s suggested additional safeguards, and this is the only simple method that is useful to find a multi-dimensional solution, which we do not require.All of the previously mentioned methods are only applicable to real data.The Lehmer-Schur algorithm is one of a number of complicated methods that is capable of solving in the complex plane, and then only for well defined polynomials (Acton, 1970).Unfortunately the Bussian equation cannot be cast in that form.Other less efficient methods include the secant method, the false position method, reiterated bracketing, Van Wijngaarden's method and Dekker's method.The Jenkins-Traub method has become fairly standard in commercially available solvers, but is extremely complicated to implement.A description of all of these methods with references is available in Press et al. (2007).
Returning to the bisection method.It is a very robust method with a long pedigree.Press et al. (2007) state that it cannot fail in the sense that it will always find one root of a single or multi-rooted function, and where there are no roots it will converge on a singularity.It converges "linearly" to a solution in the terminology of root finding algorithms, which means that convergence is mathematically exponential, or in other words, successive significant figures in the solution are won linearly with computational effort (Press et al., 2007).This rate of convergence is not bad, but may be considered slow compared to methods that converge superlinearly (i.e., improving the precision of the solution by more than an order of magnitude for each iteration).However, it does not readily overshoot a solution, which is a problem to which superlinear algorithms are prone, and which should be avoided when dealing with non-linear equations such as the Bussian equation.The concept of the method is simple.Over some interval the function is known to pass through zero because it changes sign.The method evaluates the function at its midpoint, examines its sign and replaces whichever of the intervals limits has the same sign with the midpoint.Hence the interval decreases by a factor of two for each iteration.We have used the classical bisection method as a reference with which to compare our new approach.

New approach
The Bussian equation can be written in terms of complex conductivities (σ * = σ + iσ ), complex dielectric permittivities (ε * ≡ ε − iε ), or complex relative permittivities (κ * = κ − iκ ), as in Eqs. ( 2) to (4), respectively.This is a consequence of the interchangeability of the complex conductivity and the dielectric permittivity through the relationships σ * = σ + ωε + iωε and ε * = ε − i ε + σ ω , and κ * = κ − i κ + σ ωε o , the latter of which arises due to the definition of the relative permittivity as κ ≡ εε o , where ε o is the electric constant.Note that the term ωε is equivalent to a conductivity and represents the contribution to energy dissipation made by the displacement currents, while σ represents the contribution to energy dissipation made by the conduction currents (Guéguen and Palciauskas, 1994).
The Bussian equation is non-linear as we have already discussed.The equation cannot be solved algebraically and one must use numerical methods.There is a further difficulty if the equation is to be solved for complex parameters and if m is not an integer because it then has an infinite number of roots and the trick is to be able to recognise which root corresponds to the physical solution.
Our solution of the final equation is implementable in the mathematical manipulation software Maple TM .It should also be possible to implement the method in Mathematica ® and Matlab ® , but we have not confirmed this.We have transformed the equation into a numerically solvable form using a complex number conformal mapping. where; respectively, that arise from the conformal mapping described by f : f f , and is the complex s.An analytic function is conformal mapping of a e, conformal mapping is complex variables to be and imaginary parts of it allows complex variables to be is the complex space.A conformal mapping is a transformation w = f (z) that preserves local angles.An analytic function is conformal at any point where it has a nonzero derivative.Conversely, any conformal mapping of a complex variable which has continuous partial derivatives is analytic.Hence, conformal mapping is extremely important in many areas of physics and engineering as it allows complex variables to be converted into an analytically solvable form.By letting w ≡ f (z), the real and imaginary parts of w(z) must satisfy the Cauchy-Riemann equations and Laplace's equation, so they automatically provide a scalar potential and a so-called stream function.If a physical problem can be found for which the solution is valid, we obtain a solution, which may have been very difficult to obtain directly, by working backwards.
The implementation of this solution in Maple TM code is given as Fig. 1.It is worth noting that there are only three active lines in the code, which is extremely efficient.

Testing
A number of validation tests have been carried out on the method.In all tests we have used conductivities, but note that replacement of the conductivity parameters with relative or absolute permittivities is equally valid.Hence the following results apply equally to using the Bussian equation and the conformal mapping method with relative or absolute permittivities.
In the first test we have examined that the method provides the same results as the Bussian equation for a set of limiting values.The limits of the Bussian equation by mathematical analysis are set out in Table 1, together with the results of the bisection method and the conformal mapping approach, as well as the physical limits that one might expect for a saturated rock.It should be noted that both the bisection method and the conformal mapping method provide the same limiting results as the mathematical analysis indicates except for σ m → ∞, where both the bisection method and the conformal mapping approach has difficulty in providing a root.This is more likely to be a problem with the internal solver routines of Maple TM 11 than the conformal mapping method itself, but does not impose a significant restriction because applications in this limit are negligible.The limiting solutions indicate that the results of the conformal mapping method are consistent with the Bussian equation.
The question arises whether the Bussian equation accurately describes the physical situation.The two cases where σ m → 0 and σ f → 0 are interesting because although both the conformal mapping method and classical bisection method provide the result one would expect from taking the www.solid-earth.net/1/85/2010/Solid Earth, 1, 85-91, 2010   mathematical limit of the Bussian equation, that limit is not reasonable from a physical point of view.These two limiting cases show a failing in the Bussian equation that makes it invalid at low fluid conductivities or low matrix conductivities.The failure, which is due to its derivation from effective medium theory, can be seen clearly in Fig. 2.This figure compares that Bussian equation with Archie's law (Archie, 1942), the modified Archie's law (Glover et al., 2000a) and the Korvin and Tenchov method (Korvin, 1982;Tenchov, 1998).As σ f → 0, the effective conductivity should tend towards a value defined by the porosity, cementation exponent and the conductivity of the matrix σ eff → σ m (1 − φ) m , whereas it is clear in the figure that the Bussian solution tends towards zero like the classical Archie's law.The other two curves in Fig. 2 are the solutions by Korvin and Tenchov method (Korvin, 1982;Tenchov, 1998) and by the modified Archie's law (Glover et al., 2000a).Both of these methods work better than the Bussian equation in the limit σ f → 0, while the modified Archie's law provides the exact limiting value.Note that all three models are all in fairly good agreement for σ f ≥ σ m .
In the second test we have compared the method with the classical bisection method for a loop that requires the calculation of 1001 datapoints, where a datapoint represents a set of solution parameters (i.e., [ε m , ε f , φ and m], [κ m , κ f , φ and m] or [σ m , σ f , φ and m]).In our case we chose to operate with conductivities rather than permittivities.We kept the porosity (φ = 0.2), cementation exponent, (m = 1, 1.5, 2, Figure 3. Effective electrical conductivity σ eff calculated using the Bussian equation as a function of porosity at 1001 points for both the bisection method and the conformal mapping method (φ = 0.2; m = 1, 1.5, 2, 2.5, 3; σ m = 10 -3 S/m; 10 -5 S/m ≤σ f ≤ 1 S/m).Note that the bisection method and the conformal mapping method are indistinguishable for all values of φ, m, σ m , σ f used: A graph of one versus the other is a 1:1 straight line.Fig. 3. Effective electrical conductivity σ eff calculated using the Bussian equation as a function of porosity at 1001 points for both the bisection method and the conformal mapping method (φ = 0.2; m = 1, 1.5, 2, 2.5, 3; σ m = 10 −3 S/m; 10 −5 S/m ≤ σ f ≤ 1 S/m).Note that the bisection method and the conformal mapping method are indistinguishable for all values of φ, m, σ m , σ f used: A graph of one versus the other is a 1:1 straight line.
2.5, 3.0) and matrix conductivity (σ m = 10 −3 S/m) constant, and made a calculation for 1001 different values of fluid conductivity from 1 ×10 −5 S/m to 1 S/m with twenty points per decade, distributed logarithmically.The results of the tests are shown in Fig. 3.It can be seen that both tests perform well insofar as they qualitatively produce the same results (i.e., their curves are indistinguishable).It is impossible to compare each method against a known solution, however, this test validates the conformal mapping method against a well known and respected method that is considered to be extremely robust.
On a more quantitative basis, correlation tests between the results obtained using the classical bisection method and the conformal mapping method that are given in Fig. 3, show them to be statistically the same with a covariance of (3.66 ± 6.19) × 10 −4 S 2 /m 2 , a value of 1 − r = (5.96± 1.14) × 10 −14 , and a value of 1 − r 2 = (11.92± 2.28) × 10 −14 , where r is the correlation coefficient and r 2 is the coefficient of determination.Please note that the values 1 − r and 1 − r 2 have been used because the correlation coefficient and coefficient of determination are, respectively, too close to unity to be written down effectively with precision.In these statistics, the values represent the mean calculations from each of the 5 tests for 5 different m values that are shown in Fig. 3, and the uncertainties represent the standard deviations calculated from those measurements.Applying Student's t test is a trivial exercise because the correlation coefficients are so close to unity.The similarity in the precision of the results from the two methods derives from setting values of precision and maximum number of iterations in the classical bisection method that correspond to those implicit in the Maple TM 11 solution code.This similarity of precision allows makes a direct comparison of the speed of the two methods meaningful.
Consequently, we have measured the time required to carry out the calculation of 1001 data points.The results depend upon the values of the input parameters.The results, which are given as a mean over 12 runs ± standard deviation, are given in Table 2.
The classical bisection method required between 220.69 and 280.36 s for real data depending on the value of the cementation exponent.It was fastest for m = 1, slightly slower for integer values of m and slowest when the cementation exponent was not an integer.The classical bisection method cannot be used to solve Bussian's equation if any of the input parameters are complex.By comparison, the conformal mapping method required between 3.46 and 5.41 s for real data, depending on the value of the exponent m, and between 4.02 and 27.59 s for complex data, with the time depending once again on the value of the exponent m.Hence, on elapsed time the conformal mapping method is about 52 times faster than the classical bisection method taking the most general case where m is not an integer.
Since there is a small but finite computing time associated with the structure of the program (defining variables, reporting the data etc) we have calculated the mean speed of calculation in the following way.We have run the programme for 1001 datapoints and for 1 datapoint, while retaining the same input parameters, and noting the elapsed time in each case.The mean speed was calculated by subtracting the elapsed time tor 1 datapoint from that for 1001 datapoints and dividing by 1000 (i.e., the number of datapoints calculated within   the elapsed time difference).The mean speed for the classical bisection method was 3.6 ± 0.3 and 4.5 ± 0.3 solutions per second compared to between 176 ± 16 and 289 ± 17 solutions per second needed by the conformal mapping method for real data, and between 36 ± 2 and 248 ± 24 for complex data.Hence for real data the conformal mapping method is 51 times faster than the classical bisection method, again taking the most general case where m is not an integer.
It is worth noting in Figs. 2 and 3 that all methods except Archie's law converge at σ eff = σ f = 10 −3 S/m.The Bussian equation provides a reasonable solution in the range σ f ≥ σ m .However, it should be noted that in this high salinity range the results provided by the Bussian, Korvin and Tenchov and the Modified Archie's method provide different results.
Finally, Fig. 4 shows an example of the conformal mapping method being used to solve the Bussian equation with complex input parameters.In this figure the matrix conductivity is complex σ m = 10 −3 + i10 −3 S/m, the pore fluid varies in the range 10 −5 S/m ≤ σ f ≤ 1 S/m, for a porosity φ = 0.2, and for five values of the cementation exponent, m = 1, 1.2, 1.5, 2, 2.5.It is no longer true that σ eff = σ f = 10 −3 S/m.It is also worth noting that the solution of Bussian's equation in complex space suffers from the same difficulties as in the real space when it comes to fluids with low conductivities (i.e., as σ f → 0).The innaccuracies in this limit exhibit themselves both in the real and imaginary parts of the solution and become more pronounced as σ f → 0 and as the cementation exponent increases.The problem, as described previously, resides in the formulation of the Bussian equation rather than the solution method, and arises due to assumptions that are inherent in the effective medium approach.
It should be noted that the classical bisection method is fairly straightforward to implement under any of the common programming languages, but requires a greater level of programming skill than the conformal mapping implementation in Maple TM 11.
It should also be noted that all tests were carried out on a standard desktop PC (Intel Core 2 Quad 2.4 GHz, 4 core, 3 Gb RAM, Microsoft windows XP Professional) running Maple TM 11 for the conformal mapping method.The classical bisection method was implemented using the simplified code available in Numerical Recipes in C++: The Art of Scientific Computing, 3rd edition (Press et al., 2007) and Borland C++ 5.5.In both cases the tests were run after a complete reboot of the PC and with no background tasks active.Elapsed timing was carried out using timestamps for the classical bisection method and using the Maple TM 11 native timer.There were no significant usages of physical memory for either method.

Conclusions
A new, simple and elegant method for the solution of Bussian equation for the complex effective conductivity, complex dielectric permittivity and complex relative permittivity of two phase mixtures such as water saturated rocks has been developed.The implementation of this method in Maple TM 11 allows effective and swift solution of these equations.Comparison of the conformal mapping method with the classical bisection method on the same computer shows the new method to be as precise, easier to implement and about 52 times faster to run.The new method is almost as efficient when solving the equation with parameters that are complex numbers, which is something that the bisection method, and the great majority of root finding methods is not capable of doing.
verse function of f , and is the complex rves local angles.An analytic function is onversely, any conformal mapping of a analytic.Hence, conformal mapping is ng as it allows complex variables to be ( ) f z , the real and imaginary parts of is the principal branch of ζ → (ζ − 1)/ζ (1/m) .The symbol f −1 represents the inverse function of f , and cal angles.An analytic function is ely, any conformal mapping of a tic.Hence, conformal mapping is

Figure 1 .Fig. 1 .
Figure 1.Code required to run the conformal mapping solu

Table 1 .
Limiting values of the Bussian equation solutions.