Effective transport coefficients of anisotropic disordered materials

A novel effective medium theory for homogenized transport coefficients of anisotropic mixtures of possibly anisotropic materials is developed. Existing theories for isotropic systems cannot be easily extended, because that would require geometric characterizations of anisotropic connectivity. In this work anisotropic connectivity is characterized by introducing a tensor that is constructed from a histogram of local percolating directions. The construction is inspired by local porosity theory. A large number of known and unknown generalized effective medium approximations for anisotropic media are obtained as limiting special cases from the new theory. Among these limiting cases the limit of strong cylindrical anisotropy is of particular interest. The parameter space of the generalized theory is explored, and the advanced results are applied to experiment.

Much attention was paid in recent years [1,11] to the study of certain additive quantities [27] known as intrinsic volumes, quermaßintegrals [28] (p.240), valuations [29] (p.173) or Minkowski functionals [30,31]. Additivity means that adding the values for two samples equals the sum of the value for the union plus that for the intersection [27], Eq. (42). In particular the additive Euler characteristic [27] (p.31), [28] (p.237) is often believed to determine transport in disordered media [32,33] or to "account for connectivity in a macroscale way" [1](p.325). On the other hand it is widely appreciated, that transport in random systems depends upon the existence of percolating paths from inlet to outlet [34][35][36][37][38]. Rigorous relations between the Euler characteristic and percolation were established in Eqs. (48)- (52) in Ref. [39], p. 386. Except for certain special geometries these relations show that in general the Euler characteristic cannot possibly suffice to characterize connectivity a e-mail: R.Hilfer@icp.uni-stuttgart.de (corresponding author) or transport. Most importantly, the Euler characteristic is additive, while connectivity, connectedness or "percolativity" are not.
Geometric quantities for anisotropic media, especially vectorial [53] and tensorial [54,55] extensions of additive functionals are currently under active investigation [56]. Let us emphasize that higher order Minkowski functionals as well as their tensorial extensions are studied primarily as descriptors of geometry and have, to the best of our knowledge, not been used as predictors for physical transport. Our intention in this work is to introduce novel effective medium approximations to predict percolation and conduction in disordered materials. Recall that classical effective medium theories [35,[57][58][59][60] contain only volume, more precisely volume density i.e. volume fraction, as a geometrical descriptor. In self-consistent effective medium theories for mixtures with infinite material contrast a percolation transition emerges at volume fraction 1/3. Anisotropic variants of effective medium theory based on volume fraction [41,61,62] exhibit the same percolation transition and they continue to attract much interest [42,52,63]. Minkowski functionals other than volume, however, have, as far as we know, not been used to generalize effective medium theories to anisotropic media, and the same applies to pair connectivity [38], the pair connectedness [64], or n-connectedness [65].
The self-consistent effective medium approximation for the effective dc-conductivity of homogeneous and isotropic media with porosity φ reads [35,57,58] φ (σ p − σ ) σ p + 2σ where σ p and σ m are the dc-conductivities of two isotropic materials filling pores p and matrix m (they will be introduced more precisely in Eq. (7) below). For infinite contrast σ p /σ m = ∞ this equation has a singularity at φ = φ c = 1/3 which marks the percolation transition. Nonzero solutions for σ exist only for φ > φ c , while for φ ≤ 1/3 the effective conductivity σ = 0 vanishes. Equation (1) has been generalized to anisotropic media in the special case where p is a homogeneous mixture of aligned oblate ellipsoids (x/b) 2 + (y/b) 2 + (z/c) 2 ≤ 1 with half axes b ≥ c > 0. In this case the effective conductivity σ becomes tensorial. Its eigenvalues σ x , σ y , σ z obey σ x = σ y and Eq. (1) becomes essentially two coupled equations [41], Eq. (8) for σ x = σ y and σ z . The depolarization factors N a , N c depend on a combined ratio R of the two unknowns σ x , σ z and the two half axes b, c through the relations that couple the two Eqs. (2a) and (2b). As in the isotropic case, these two coupled equations exhibit a percolation singularity at φ = φ c = 1/3. Local porosity theory [36] is a generalized selfconsistent effective medium theory that contains local porosity distributions and local percolation probabilities above and beyond volume fractions to characterize the geometry and connectedness of a disordered medium. An important consequence are percolation thresholds φ c = 1/3. Predictions from local porosity theory agree with observations such as the empirically known Archie's law [36]. Local porosity distributions and local percolation probabilities have become standard tools for image analysis of porous media [66,67]. They were applied to clay rocks [68], fractured coal [69], Opalinus clay [70,71], chalk [72], reticulate porous ceramics [73], macroporous reticulated silicon oxycarbide [74], dielectric characterization of wood [21], highly porous polymeric matrices [75], metallic filters [44], beryllium pebbles and arctic firns [64], and recently even to Ibuprofen tablets in pharmaceutical research [76]. Numerous publications [17,32,[77][78][79][80][81] have demonstrated the predictivity of percolation probabilities and the total fraction of percolating cells for physical transport.

Problem and objective
A shortcoming of local porosity theory for electrical transport in heterogeneous media (introduced in Ref. [36] and reviewed in Refs. [37,82]) is that the local percolation probability λ, which characterizes the geometrical connectivity of the medium, is a scalar but not a tensorial quantity. Disordered or heterogeneous media, however, generally have tensorial transport properties even if the constituent materials are isotropic.
An attempt to extend LPT to anisotropic systems was made in Ref. [83]. The anisotropy is described using eigenvalues of the "orientation matrix". The orientation matrix is constructed from the position vectors of those voxels within a measurement cell, that belong to the largest percolating cluster in the cell. The position vectors are computed relative to the "cluster center" defined as an arithmetic average. Two eigenvalue ratios of the orientation matrix are collected into local anisotropy distribution. The relation of the local anisotropy distribution to local percolation probabilities is not discussed in Ref. [83]. As a consequence the generalized effective medium equation is unchanged and transport quantities such as conductivities or permeabilities remain scalar quantities.
The objective of this work is to extend local porosity theory and effective medium theories to anisotropic media in such a way, that the transport quantities become tensorial. For brevity and simplicity we consider only transport of charge, i.e. the case of electrical transport.

Model
Let p, m ⊂ R d denote the domains (pore space, matrix) occupied by the two components of a heterogeneous or porous sample S = p∪m ⊂ R d . The Maxwell equations [84] read for r ∈ p ∪ m, where σ is the electrical conductivity tensor. The symbols 0 , r denote the of dielectric permittivity of the vacuum resp. the permittivity tensor of the medium, and μ 0 , μ r are the scalar resp. tensorial magnetic permeabilities of the vacuum resp. the medium with μ 0 = 4π × 10 −7 N/A 2 and 0 = 1/(μ 0 c 2 ) ≈ 8.8542 × 10 −12 F/m. The conduction current due to free charges J c appearing in Ohm's law, Eq. (4a), is related to the total current via the current of bound charges defined as J b (r, t) = J(r, t) − J c (r, t). Fourier transformation with respect to t gives where ω = 2πν is the angular frequency and ν denotes frequency. One has (r, ω) = (r, ω) + i (r, ω) for the complex conductivity and dielectric tensors. The material parameter functions for a composite material consisting of two homogeneous but anisotropic components are assumed to have the form is the characteristic (or indicator) function of a subset X ⊂ R d . The relaxation frequencies, defined as are assumed to exist. They can be used to make ω dimensionless.
The discontinuity of the material parameters at the interface requires one to specify boundary conditions at ∂p = ∂m. Mathematically the equations are interpreted in a weak sense as equations for distributions. Depending on the boundary conditions a suitable domain could be a Sobolev space for the potentials resp. a space of potential fields for the electromagnetic fields [85,86]. The boundary conditions at the interface are for r ∈ ∂p = ∂m, where Q ∂ (r, ω) (resp. J ∂ (r, ω)) are the Fourier transforms of (possibly time dependent) surface charge (resp. surface current) densities with support in ∂p = ∂m. The notation B ∂p (r, ω) (resp. B ∂m (r, ω)) is the limiting value of the vector field as the point r ∈ ∂p = ∂m is approached from within p (resp. from within m). It is assumed that there are neither volume nor surface charges or currents inside the medium so that Q = 0, Q ∂ = 0, J = 0 and J ∂ = 0 from now on. This assumption is commonly made to exclude phenomena and effects arising from electrical double layers and other physico-chemical processes in the interfacial region.
When the wavelength is large compared to the scale of heterogeneities the time derivatives in in Faraday's law, Eq. (3), are small compared to spatial derivatives. Setting them to zero results in the quasistatic approximation valid for small frequencies where denotes the scale of heterogeneities and c is the speed of light. For heterogeneities with ≈ 10...100μm one finds ω 3...30THz approaching the infrared. The quasistatic approximation decouples the full equations (3). Fourier transformation of eqs. (3) with respect to t results in the quasistatic equations where r ∈ p ∪ m in Eqs. (12a)-(12c), r ∈ ∂p = ∂m in Eqs. (12d)-(12e), and where (r, ω) denotes the complex frequency dependent local dielectric function. The effective macroscopic dielectric tensor of selfconsistent effective medium approximations [58] is usually defined by averaging Eq. (12c) where the angular brackets · denote (spatial or ensemble) averaging and Π are parameters of physical importance (such as porosity or connectivity) that arise from averaging the smaller scale heterogeneities.

Disordered media
A disordered (or porous) medium m (or p) is assumed to be a random set. A random set is a set-valued random variable [87]. The probability space of random disordered geometries is denoted (Ω, A, P ). The set Ω is the set of all closed subsets of R d including ∅. The σalgebra is generated by the system of closed sets having nonempty intersection with compact sets and the random geometry distribution P is the image measure of the measure on the probability space underlying the setvalued random variable. For mathematical details on random sets the interested reader is referred to Ref. [82].

Local porosity
Let K ⊂ R d be a convex and compact set with centroid at the origin 0 ∈ R d . The centroid is the geometrical center of K. Then denotes its translate by a vector r ∈ R d . The local porosity in a measurement cell K(r) placed at position r is defined as where is the volume of a set K. The local matrix volume fraction is φ m = 1 − φ.

Anisotropic local matrix fraction
Consider from now on d = 3 and let m i denote the local matrix space clusters defined as (path-)connectedness components of m ∩ K(r). Then where m i ∩ m j = ∅ for i = j and n m (K(r)) is the total number of local clusters within K(r). The barycenters of these local clusters are and is their local volume fraction. The matrix is assumed to have a constant density taken to be unity. Averaging over the local clusters the tensor is introduced, where X(r) = r ⊗ r, is the tensor product. The tensor is real, symmetric and non-negative definite. In the following the tensor is normalized as

Local percolation
Let ∂K(r) denote the boundary of a measurement cell with centroid r. A measurement cell K(r) = r + K is called percolating if there exists a continuous path with p 0 ∈ ∂K(r) (21b) Fig. 1 Schematic illustration from Ref. [88] of the local percolation criterion (21). The grey hatched circle centered at p 0 is the forbidden ball B (p 0 , |r − p 0 |). Two paths are shown that start at p 0 on the left boundary. The dashed path is percolating according to the percolation criterion (21), because it ends outside the forbidden ball. The dashdotted path is not percolating, because it ends inside the forbidden ball where B(x, a) = {r ∈ R 3 : |r − x| ≤ a} denotes a ball of radius a centered at x. Equation (21) is illustrated in Figure 1. The set of local percolating directions for a measurement cell K(r) centered at r is defined as where S 2 = ∂B(0, 1) is the unit sphere and p 0 p 1 means that p 0 is the starting point and p 1 is the end point of a percolation path in the sense of Eq. (21). The local percolation indicators, defined as are a family of {0, 1}-valued random variables indexed by K. The image of the measure P under the map defines the total percolation probability p(K) of the measurement cell K.

Anisotropic local percolation
The set L of local percolating directions defines a tensor with respect to the cell center r. It is again normalized by dividing with the trace. The non-negative eigenvalues J pa (K(r)), J pb (K(r)), J pc (K(r)) obey J pa (K(r)) +J pb (K(r)) + J pc (K(r)) = 1, and are assumed to be ordered as 1 ≥ J pa ≥ J pb ≥ J pc ≥ 0.

Local geometry distribution
The local geometry of a random set p inside a measurement cell K is characterized by the geometric parameters g(p; and parametrize the anisotropy of the measurement cell K with local porosity φ. These quantities form a family of G-valued random variables of P under g is the family of local geometry distributions with family index K.
In practical applications the distribution P of the random geometries is unknown. Instead one is given a sample assumed to be a realization of a random geometry. The local geometry distribution is then estimated as the joint empirical probability measure where r i is a sequence of cell centers such that all K(r i ) are mutually disjoint, and N denotes the total number of different placements of K. Mathematically this empirical density is an estimator for a Young measure (see Refs. [39,82,89] for details).

Limiting local geometry distributions
The dependence of μ(g; K) on K disappears in the limit where K shrinks to a point or expands to become R 3 . In the limit |K| → 0 of small pointlike measurement cells one has with Eq. (27) where φ = φ(S) is the total porosity of the sample. In the opposite |K| → ∞ limit the limiting local geometry distributions are also expected to become independent [89] of K. Then Eq. (27) gives where μ 1 (φ, d p , e p ) is the conditional probability density that a cell is percolating in the scaling limit, has porosity φ and percolation anisotropy parameters d p , e p , while μ 0 (φ, d m , e m ) is the conditional probability density that a cell is non-percolating in the scaling limit, has porosity φ and percolation anisotropy parameters d m , e m . The quantity λ(φ) is the conditional probability that a very large measurement cell is percolating given that it has local porosity φ. It may be called limiting local percolation probability (cf. Ref. [36]) in the sense is expected to hold from Eq. (25). The possible macroscopic limits have been classified in Ref. [89], Sec.VII.b, into four cases. In the following the macroscopically heterogeneous cases will be neglected and the porous medium is assumed to be macroscopically homogeneous and anisotropic. If it is macroscopically non-random one has where φ is the macroscopic porosity and d p , e p , d m , e m are the macroscopic anisotropy parameters for pore and matrix. The simple choice will be used below. Expectation values of measurable functions f : G → R are computed as

Formulation of the local problem
Let d = 3 and let r T = (x, y, z) ∈ R 3 denote a vector in cartesian coordinates. The potential problem for the unknown electrostatic potentials U c , U s , U u in the domains c, s, u ⊂ R 3 reads where c , s , u are the dielectric tensors of the homogeneous but anisotropic dielectric materials filling the domains c, s, u and E ext denotes a constant applied field. The union of c and s forms an ellipsoidal inclusion consisting of an ellipsoidal core surrounded by a shell that is embedded into an environment (or background) defined as the complement of e. It is assumed that the quadratic forms Q e , Q c are non-degenerate, and that c ⊂ s and s \ c = ∅ holds throughout. The finite twodimensional boundaries are denoted and the vector fields n ∂e , n ∂c in Eq. (34) are unit normal vector fields on e and c. The boundary conditions (34d), (34e) at ∂u∩∂s = ∂e and (34f), (34g) at ∂s ∩ ∂c = ∂c express continuity of the potentials and the normal component of the dielectric displacement. The boundary condition (34h) demands that the potential vanishes at infinity |r| → ∞. The non-degenerate quadratic form matrices Q e , Q c (with det Q e = 0, det Q c = 0) depend on the material properties s of the shell as follows. Assume that det s = 0 and define the coordinate transformation √ s by the formula as the square root of s . The linear, non-orthogonal and non-singular coordinate transformation (36b) relates the vector r = (x, y, z) T in the original coordinate system with the vector r s = (x s , y s , z s ) T in the transformed coordinates. To specify the transformation uniquely it is assumed that holds. Then a c,s < a e,s , b c,s < b e,s and c c,s < c e,s and the coated ellipsoid is parameterized by four parameters. If the half axes a e,s , b e,s , c e,s are normalized by division with their sum a e,s + b e,s + c e,s , then three parameters suffice to characterize the coated ellipsoid. They could be chosen as the volume fraction φ c , and two normalized half axes, or as φ c and two of the axis ratios for shell and core. Here the axis ratio ϕ c is a core to shell length fraction for the intermediate axes. The confocality condition (38) implies the relations between these axis ratios. The core c within the coated ellipsoidal inclusion e has a local volume fraction that can be used to eliminate ϕ c in exchange for φ c .
Given d s , e s and φ c these equations are solved for ϕ c to determine as functions of d s , e s and φ c . The functions D c , E c : 1] show clearly, that d c , e c are not independent parameters. Instead, they are fixed as soon as φ c , d s , e s are given.

Solution of the local problem
To solve Eq. (34), one first transforms Eq. (34c) into the Laplace equation in vacuum using the transformation where it is assumed that det u = 0. To specify the linear, non-orthogonal and non-singular coordinate transformation √ u uniquely it is assumed that Following Refs. [90] and [91], Eq. (44), the solution of the boundary value problem (34) can be written as where ξ u = ξ u (r) (resp. ξ s = ξ s (r)) is the first ellipsoidal coordinate of the vector r in (x u , y u , z u )-coordinates (resp. (x s , y s , z s )-coordinates) and C > 0 is the confocal constant in Eq. (38). The ellipsoidal coordinates Similarly, the ellipsoidal coordinates (ξ s , η s , ζ s ) of a point (x s , y s , z s ) are defined as the three roots of the equation where |c| is the volume of c and w = a c,s , b c,s , c c,s for i = x su , y su , z su . The diagonal entries N where w = a e,s , b e,s , c e,s for i = x su , y su , z su and |e s | = |e| is the volume of e. Then are the depolarization tensors in (x, y, z)-coordinates.
With the abbreviations N e,u = N e,u (0), L e,u = L e,u (0), where φ c = |c|/|e| is the local core volume fraction and is the dielectric difference between the core and the shell material.

Self-consistent approximation
The potentials U c (r; c , s , u ), U s (r; c , s , u ), and U u (r; c , s , u ) in Eq. (45) depend on the data c , s , u . The dielectric permittivity e of the coated ellipsoidal inclusion is determined as the solution of the equation for all r ∈ u by demanding that the potential outside of the inclusion is identical to that of the full solution.
Setting c = s = e in Eq. (52) gives for the potential and for the fields. Equating Eq. (54b) with Eq. (52c) yields where and Solving Eq. (56) for e gives which simplifies to become for the dielectric permittivity of the ellipsoidal inclusion in terms of its geometric and material parameters. The geometric parameters are restricted to φ c , d s , e s because L c,s follows from L e,s and the latter depends only on d s , e s . Following Ref. [58] and Eq. (13), the effective dielectric permittivity is defined by averaging Eq. (55) over the disorder to get where . . . is the disorder average. The implicit dependencies of e are clear from Eq. (60). The implicit dependencies of L e,u emerge from noting that the quadratic form matrix Q e of the ellipsoidal inclusion from Eq. (35a) can be expressed equivalently as in r u -coordinates or r s -coordinates. Noting that Eq. (48) can be written as du (63) allows to rewrite L e,u as ds (64) and this shows that L e,u = L e,u ( u , s ) depends on s and u . Self-consistency determines and L such that the background equals the effective medium. Setting the self-consistent effective medium approximation for . The effective depolarisation tensor L defined as with Q −1 e given in Eq. (62b) depends on . In addition to also the shape of the ellipsoidal inclusion e has to be determined self-consistently. This differs from conventional Bruggemann theory.

Local porosity theory
The self-consistent equation of local porosity theory is obtained by expressing e and L in terms of local geometric quantities and replacing the disorder average with an average over the local geometry distribution μ. This gives with . . . μ given in Eq. (33). The dependencies have to be related to the materials filling p and m.
The cases of Λ = 0 (nonpercolating local geometry) and Λ = 1 (percolating local geometry) need to be distinguished, because for percolating local geometries the shell s should contain p-material, while it should contain m-material in the nonpercolating case. The core and shell properties are specified as Inserting Eq. (69) into Eq. (68) and using (33) gives (70) as the mixing law of local porosity theory. The unknown quantity in Eq. (70) is . The functional e is given in Eq. (60) with Eqs. (57) and (52d), while L is from eq.(67).

Small measurement cells
In the limit |K| → 0 of pointlike measurement cells the mixing law of local porosity theory becomes after inserting Eq. (30) into Eq. (70). For isotropic materials p = p 1, m = m 1 the classic Bruggemann effective medium approximation (EMA) equation ensues, which should be compared to Eq. (1). The known anisotropic effective medium formula (2) is recov- where b, c are the parameters in Eq. (2). The limit of small measurement cells is the limit in which the basic assumption of statistical independence underlying effective medium approximations is violated. It is also worth to remark that, in a limiting sense, Eq. (2) requires measurement cells to remain anisotropic when they shrink to a point.

Large measurement cells
For large measurement cells |K| → ∞ Eq. (32) gives For geometrical isotropy N c,s = N e,s = 1/3 and material isotropy p = p 1, m = m 1 this becomes the generalized effective medium approximation from Ref. [36], Sec.V.A, where are the local dielectric functions.

Material isotropy
A situation of great practical and applied interest is the case of material isotropy but geometrical anisotropy. In this case the effective material parameters the effective dielectric properties become anisotropic only due to geometric properties of the mixture, while the materials themselves are isotropic. General results for macroscopically random systems will be investigated in the first Sect. 9.1 of this section. The non-random case will be discussed in more detail in the remaining subsections.

Macroscopically random systems
Isotropic materials occupying p and m have dielectric tensors proportional to 1 where the scalar prefactors are in general frequency dependent. Then c = c 1, s = s 1 so that where N e,s , N c,s and hence L e,s , L c,s are all diagonal. It follows now from Eq. (62b) that Q −1 e is diagonal and hence from Eq. (62a) that u is diagonal, and thus that L e,u is diagonal. As a result and L are diagonal. Equation (60) becomes where the diagonal matrices N c,s , N e,s , given in Eqs. (49) and (50), can be written as  (67) gives

.(90c)
This form is convenient, because it allows one to read off the limiting cases for all 0 ≤ d, e ≤ 1.
Compared to the elliptic integrals in Eqs. (48), (49), (50) the parameters d s , e s in Eq. (88b) are neither from the unit interval nor real, but in general complex numbers. As a consequence the well known formulae for real numbers given in Ref. [93] do not apply. The results for N a , N b , N c with complex parameters d, e have been computed below in Appendix A.

Macroscopically non-random systems
For macroscopically homogeneous but non-random systems Eq. (32) yields the diagonal equation is used again for local functions. Written in components the diagonal Eq. (92) becomes three coupled equations (j = x, y, z) where the abbreviations are the diagonal elements in N e,s , N c,s from Eq. (84), and p , m are from Eq. (80). Here from Eq. (42). The material parameters p , m , σ p , σ m in p , m and the geometrical parameters φ, λ, d p , d m , e p , e m in these equations are measurable or known from experiment. The complex effective depolarization factors for p and m depend on . They couple the three equations (94) via Eq. (88b). Clearly, j = 0 solves Eq. (94) for any j. For j = 0 one can rewrite the system of equations as where the normalized effective dielectric permittivity is the unknown, is the right hand side, and with j = x, y, z are parameters. Equation (97) can also be formulated as This is not a solution of Eqs. (94) but a reformulation, because A j , B j depend through R j on the two ratios x / y and y / z . To obtain analytical information about solutions it is useful to consider limiting and special cases. The identity j=x,y,z holds generally. It is obtained by summation of Eq. (97a) with the help of Eq. (A13).

Spherical (Isotropic) Limit
The spherically isotropic limit is the limit d p → 1, and d m → 1, and e p → 1, and e m → 1. In this limit N e,s = N c,s = 1/3 and N pj = N mj . Inserting this into Eq. (97) renders P j , M j independent of j and R j = −N pj follows from Eq. (95a) and (95b). Subtracting Eq. (97a) with j = x from (97a) with j = y gives where Eq. (97c) was used. The integral on the right hand side depends on z , while the left hand side does not. Thus the equation can only be fulfilled for x = y . Repeating the same subtraction with j = y and j = z gives y = z , now invoking that the integral depends on x . In this way the equality and the generalized effective medium equation from Ref. [36] is recovered.

Planar (flat) limits
The  N C x , N C y , N B x , N With this (97d) and (97e) give and eqs. (97a) simplify to become for j = x, y, z. Note that these equations do not decouple. For d p = d m and e p = e m Eq. (103) is independent of λ and hence the solution becomes independent of λ.

Percolation limit
For sufficiently small frequencies ω ω p , ω ω m the effective dielectric function j (ω) → σ (0) approaches the effective d.c. conductivity and the same holds for the material parameter functions p (ω) → σ p (0) and m (ω) → σ m (0). The complex quantity K approaches a non-negative real value which will be called conductivity contrast. The d.c. conductivities obey σ p (0), σ m (0) ≥ 0. Without loss of generality it can be assumed that holds true, i.e. that the material of higher conductivity occupies the pore space p. Vanishing contrast κ = 0 means σ p = σ m . In the limit κ → ∞ of infinite conductivity contrast one has for j = x, y, z in Eq. (97d) and (97e). Then Eq. (99) becomes is the normalized d.c.-conductivity. Because σ j ≤ P ∞ j , the left hand side is non-negative. For λ ≤ 1/3 the right hand side becomes negative, and hence λ = 1/3 is a percolation threshold.
Because M j = 0, Eq. (98c) gives B j = 0. Thus σ j = 2A j and in the percolation limit. These are three coupled equations, that can be reduced to two coupled equations by taking ratios. They have to be solved for the two unknown ratios u xy , u yz . The functions R x , R y , R z depend only on these two ratios. For λ = 1 and for λ = 0 one finds because the expression P ∞ j N pj /(N pj − 1) resulting from Eq. (109) becomes negative.

Planar (flat) LPT
In the planar case with e p = e m = e → 0 the ellipsoids degenerate into flat ellipses. As discussed in Section 9.4 one has d for the local functions. Then Eq. (109) implies that σ z = 0 for all λ. Furthermore, σ x = σ y = φ must hold for λ = 1, because of Eq. (111).

Ellipsocylindrical LPT
In this case d = 0 and e p = e m = e and 0 < e < 1.
Evaluating the integrals in Eq. (90) for d = 0 gives eqs. (91f) resp. (C5). From Eq. (40) one finds so that Eqs. (106) and (96) yield for the local functions, and holds true for their ratio. Dividing Eq. (109) for j = y by Eq. (109) for j = z results in the simple expression for the ratio of the unknown effective conductivity components. Inserting this into Eq. (95a) gives the surprising result for all e, i.e. independent of σ j or e. The result decouples eqs. (109) to give with the help of Eq. (106) Eq. (127). For λ ≤ 1/2 it follows that σ y = σ z = 0 and inserting this into Eq. (107) gives

Percolation transitions
The Bruggemann effective medium approximation (1) exhibits a percolation transition at porosity φ c = 1/3 at which the dc-conductivity σ (0) vanishes. While qualitatively correct, the percolation threshold is generally too high for direct application to porous and disordered media. They often become percolating at much lower volume fractions. Indeed, Archie's empirical law [36], which relates electrical conductivity to porosity, suggests that there is no lower limit for the percolation threshold. The Bruggemann result (1) has been generalized from spheres to oblate ellipsoids in Ref. [41]. It is governed by the generalized Eq. (2) which reduces to Eq. (1) for b = c. The numerical solution of the anisotropic or ellipsoidal Bruggemann equation (2) is shown as dashdotted lines in Fig. 2 as a function of φ. The isotropic or spherical Bruggemann result is plotted as the solid line Fig. 2. The percolation threshold is φ c = 1/3 in both cases.
The numerical solution of Eq. (109) for a comparable oblate case with (d p , e p ) = (1.0, 0.5) is shown as two dashed lines in Fig. 2. The long dashed line is σ x = σ y and σ z is short dashed. Here, λ(φ) = φ 1/3 has been assumed for the local percolation probability, and the remaining two axis ratios are (d m , e m ) = (1.0, 0.01). The percolation threshold φ c is then found as the solution of the equation λ(φ) = 1/3, which gives φ c = 1/27 ≈ 3.7%. This is confirmed in Fig. 2. It is much more realistic than φ c = 1/3 for natural media. The value of φ c depends on the local percolation probability λ(φ), thereby highlighting the importance of this geometrical quantity. The solutions of Eq. (2) and Eq. (109) both display anisotropic σ conductivity. In both cases σ x = σ y , because the anisotropy is of oblate type. However, several differences exist between Eq. (2) and Eq. (109).
While σ x,Sch and σ z,Sch in Fig. 2 do not differ much, the difference between σ x and σ z is significant for φ < 1. This is in part due to the small value of e m . In addition σ z does not approach unity for φ → 1. This behaviour is caused by the local functions and stems from the confocality condition. It appears in the limit φ → 1 whenever e p = 1 and is absent in the prolate case when e p = 1. Of course, the two additional axis ratios (d m , e m ) in Eq. (109) allow for a greater variety of solutions than for Eq. (2). Recall also that Eq. (109) is only one of the many special cases from the generalized LPT in Eq. (70).

dc-Conductivity at infinite contrast
This section gives numerical solutions of the anisotropic effective medium equations (109) for the unknown normalized dc-conductivity σ (0). Equations (109) hold for material isotropy and geometric anisotropy in the percolation limit, i.e. at frequency ω = 0 and infinite conductivity contrast κ = ∞. The solutions for the unknown real part of the normalized effective conductivity σ are displayed in twelve subfigures of Fig. 3.They depend on six parameters (φ, λ, d p , e p , d m , e m ) with values in the six dimensional unit cube [0, 1] 6 . The porosity is fixed at φ = 0.1 in all cases.
The central subfigure in Fig. 3   The subfigures (a), (d), (g) and (j) located on the four corners of Fig. 3 illustrate cases with d p ≈ d m and e p ≈ e m . In the isotropic case (a) the three components σ x , σ y , σ z are nearly identical as expected from Eq. (101). In the circular ---circular case (d) one has σ x ≈ σ y but σ z ≈ 0. In both these cases the percolation transition occurs at λ = 1/3. In the flat --flat case (g) a second percolation transition starts to appear in σ x and σ y at λ = 1/2 while σ z ≈ 0 remains small. Note also that σ x and σ y are nearly linear above λ ≈ 1/2. In the circocylindrical ---circocylindrical case (j) the percolation transition in σ x at λ = 1/3 and at λ = 1/2 in σ y and σ z remain clearly visible. Now all components are nearly linear above λ ≈ 1/2. The main difference to the flat ---flat case (g) is that σ z = 0 and σ x = σ y at λ = 1. Other cases with d p ≈ d m and e p ≈ e m interpolate between these four limiting cases. and (k). In all these cases the percolation transition in σ x occurs at threshold 1/3, while σ y approaches a percolation transition with threshold λ = 1/2 as d p → 0.

Exchanging
The parameter pairs (d p , e p ) and (d m , e m ) have very different impact on the result. This is illustrated in an approximate way by exchanging them.
Subfigures (b) and (c) illustrate such an exchange for parameters near the prolate and oblate boundary. Note that the curvature of σ y changes.
Subfigures (e) and (f) illustrate the exchange between positions close to the flat boundary and positions somewhat removed from the upper prolate boundary. Again the curvature of σ y changes, but in addition the change in σ z is significant. In the case (f) the component σ z is nearly zero for all λ, while it becomes strongly curved upwards in subfigure (e). This indicates a percolation transition at λ = 1 in the limit e m → 0. Keep in mind, that the condition λ = 1 can be fulfilled for any 0 ≤ φ ≤ 1.
The subfigures (h) and (i) illustrate the exchange between points near the flat and cylindrical boundary. Subfigure (h) shows that σ z ≈ 0 holds for all λ. There is no percolation transition in σ y at λ = 1/2 in (h). In (i) a percolation transition in σ y occurs at λ = 1/2, while the strong curvature of σ z indicates again a percolation transition at λ = 1 for e m → 0.
The cylindrical cases with d p ≈ 0 or d m ≈ 0 are shown in subfigures (k) and (l). In (k) the percolation transition for σ y and σ z at λ = 1/2 is clearly visible. This transition does not appear in subfigure (l).
The eight subfigures illustrate the difference between the parameter pairs (d p , e p ) and (d m , e m ). It is mainly the pair (d p , e p ) that determines whether or not an additional percolation transition emerges at λ = 1/2. It is also responsible for the smallness of σ z in the limit d p → 0 or e p → 0.

Interpolation
The twelve results for σ shown in Fig. 3 depend in many cases continuously on the parameters (d p , e p , d m , e m ). The continuous parameter dependence can be used to interpolate between different cases and to envisage results for other cases. This will be illustrated with two examples.
Consider first the sequence (d) → (f) → (h) → (j). In the first transition the pair (d m , e m ) moves from the lower boundary (flat) to the upper boundary (prolate) while (d p , e p ) moves slightly away from the oblate boundary. The change in σ is small and mainly seen in σ y . In the transition (f) → (h) both parameter pairs move towards the cylindrical boundary. The result is that the curvature in σ y becomes much more pronounced and also σ x becomes curved. The curvature in σ x is due to the percolation transition that emerges at λ = 1/2, because the open symbol has moved to the left. This is confirmed in the transition (h) → (j) as the parameter pair (d p , e p ) moves even closer to the cylindrical boundary. The change in σ y on the other hand is caused by the closed symbol being close to the cylindrical boundary.
Next, consider the sequence (c) → (e) → (i) → (l). In the first transition (c) → (e) the pair (d p , e p ) remains essentially fixed, while (d m , e m ) moves from prolate to flat. As a result σ z becomes strongly curved upwards, indicating a percolation transition at λ = 1. Shifting (d p , e p ) close to the cylindrical boundary and (d m , e m ) only slightly keeps the upward curvature in σ z and exposes the emerging percolation transition at λ = 1/2 in σ x and σ y . The transition (i) → (l) then undoes most of these changes and suggests to close the circle with an additional transition (l) → (c).  94) for the normalized frequency dependent conductivity tensor σ (ω) at infinite conductivity contrast κ = ∞. The line styles for σ x , σ y and σ z in Fig. 4 are the same as in Fig. 3. Each of the subfigures in Fig. 4 has a fixed parameter quadruple (d p , e p , d m , e m ) and is labeled (a), (d), (g), (j), (e), and (i), to indicate the analogous pair (d p , e p ) ----(d m , e m ) labelled identically in the central subfigure of Fig. 3. The six subfigures in Fig. 4 are further subdivided into four subfigures. The two columns are distinguished by λ = 0.35 (left column) and λ = 0.55 (right column). The values are chosen in this way, because they are slightly above the two percolation thresholds at λ = 1/3 and λ = 1/2. In each subplot the upper two plots show log 10 ( σ ) against log 10 (ω), the lower row displays the behaviour of the derivative d log 10 ( σ (ω))/d log 10 (ω) against log 10 (ω). The numerical value of the derivatives may be viewed as a local scaling exponent. The abscissa in all plots is log 10 (ω). The frequency is dimensionless and given in units of the relaxation frequency ω p defined in Eq. (9).

ac-Conductivity at infinite contrast
The porosity is φ = 0.1 in all cases.

General observations
Two different sources can be identified which cause dispersion of σ . The first source is the dispersion of the components P j and M j , defined in Eq. (97d) and Eq. (97e). It is the dispersion in the response of a single local geometry. The second source causing dispersion arises from disorder and the percolation limit. The spherical case (a) is shown for reference in the upper left corner of Fig. 4. In this case the dispersion is weak and it has a limited range in frequency.
The ordering σ x (ω) ≥ σ y (ω) ≥ σ z (ω) holds in all subfigures as for dc-conductivities. In general, the dispersion decreases for increasing values of λ. Contrary to the limit ω = 0, the values of σ at high frequencies are barely affected by λ.

Flat pores e p ≈ 0
The case (d) from Fig. 3 is shown in the upper right corner of Fig. 4. Clearly, σ z σ x and σ z σ y holds for all frequencies. As in the case ω = 0, this results from P z → 0 for e p ≈ 0. For σ x and σ y , there is one dispersive frequency range, which arises from the percolation transition at λ = 1/3. For σ z , there is an additional dispersion at higher frequencies. It arises from the dispersion of P z .

Cylindrical pores d p ≈ 0
The case (j) from Fig. 3 is placed on the right side of the middle row in Fig. 4. As discussed in section 10.2.3, σ y and σ z approach a percolation transition at λ = 1/2. This yields two effects for the dispersion of these components at λ = 0.35. First, the dispersion is much larger than the dispersion of σ x , and second, it occurs in a much broader frequency range. Both observations can be explained by the fact that the percolation limits at λ = 1/3 and λ = 1/2 lead to dispersion in two different frequency ranges which overlap. This effect cannot be observed at λ = 0.55 > 1/2.

Completely flat case
The case (g) from Fig. 3 is shown on the middle left side of Fig. 4. The dispersion of σ x and σ y is similar to case (j), because σ x approaches the percolation limit at λ = 1/3 and σ y at λ = 1/2. The behavior of σ z can be qualitatively understood as a combination of the cylindrical (j) and the circular case (d). Compared to case (j), there is an additional dispersion at high frequencies which arises from the dispersion of P z for e p ≈ 0.

Flat matrix e m ≈ 0
The case (e) from Fig. 3 is shown in the lower left corner of Fig. 4. At first sight the dispersion looks similar to the flat case (g). However, the additional dispersion at high frequencies for σ z results from the percolation limit at λ = 1, and not from the dispersion of P z for e p ≈ 0. This can be concluded from the fact that this dispersion barely changes between λ = 0.35 and λ = 0.55. The case (i) from Fig. 3 is shown in the lower right corner of Fig. 4. Because d p ≈ 0, the conductivity eigenvalue σ y approaches the percolation limit at λ = 1/2, and σ z approaches the percolation limit at λ = 1 because e m ≈ 0 holds. This means that the three dif-  109) for σ x (φ) (solid curve), σ y (φ) (dashed curve) and σ z (φ) (dash-dotted curve) as functions of average porosity φ for anisotropy parameters dp = ep = dm = em = 0.58 and local percolation probability λ(φ) as in Fig. 2. The percolation threshold is at φ c = 1/27 ≈ 3.7% as in Figure 2. Experimental data for shaly sands, normalized as in Eq. (108), are displayed as symbols. Square symbols are the data from Table 1, while  circles are from Table 3 in Waxman and Smits [95] ferent eigenvalues of σ approach three different percolation limits. From this follows that the components are dispersive over different broad frequency ranges. At λ = 0.35, the dispersion at low frequencies results from the transition at λ = 1/3, the dispersion at intermediate frequencies from the transition at λ = 1/2, and the dispersion at high frequencies from the transition at λ = 1.

Application to experiment
A simple application of the new approximations to experiment can be made for anisotropic shaly quartz sands. It has recently been suggested by Nguyen, Vu, and Vu [94] that the electrical conductivities of shaly sands measured by Waxman and Smits [95] and Hill and Milburn [96] are anisotropic. Anisotropy ratios as large as 5 between normal and transversal conductivities are suggested in Nguyen, Vu, and Vu [94], Table 1. For details on the samples and their description the reader is referred to Waxman and Smits [95], p. 113. The electrical conductivity of the anisotropic shaly sands was measured using an impedance bridge for only one particular orientation. The data with an accuracy of 0.1 percent will be viewed here as resulting from randomly orienting the sample. Figure 5 shows the logarithm of the measured d.c. electrical conductivities normalized according to Eq. (108). Square symbols are the data reported in Table 1, while circles are from Table 3 in Waxman and Smits [95].
Detailed sample specific information for the anisotropy of pore p or matrix m in the shaly sands are not available from the publication [95]. To apply the theory nevertheless, it will be assumed for simplicity that d p = e p = d m = e m holds true. For the local percolation probability λ(φ) the function underlying Fig. 2 is assumed again, so that the percolation threshold occurs again at porosity φ c = 1/27 ≈ 3.7%. The numerical solution of Eq. (109) yields the three eigenvalues σ x (φ), σ y (φ) and σ z (φ) of the conductivity tensor as functions of average porosity. The three curves for d p = e p = d m = e m = 0.58 are displayed in Fig. 5.
The figure suggests that the scatter in the experimental data may indeed arise from an anisotropy of the samples. In conclusion, the theoretical assumptions underlying Eq. (109) and the experimental indications, that the measured conductivities of shaly sands are anisotropic, can be consistently quantified with the help of the approximations developed in this paper.

Conclusion
Generalized effective medium approximations for anisotropic media have been derived from local porosity theory. A purely geometric and experimentally measurable characterization of the anisotropy in (local) connectivities was introduced as a prerequisite for the generalized theory. The resulting anisotropic local porosity theory is based on the familiar local geometry distributions from the isotropic case. These distributions characterize the complex pore-space geometry. The generalized theory contains isotropic local porosity theory as a special case. The local geometry distributions characterizing local anisotropic connectivity are obtained from local percolating directions. Local percolating directions were defined based on the generalized percolation criterion introduced in Eq. (21). The geometric quantities enter directly into the generalized anisotropic selfconsistent effective medium equation. No adjustable parameters appear in the generalized anisotropic local porosity theory. Several hitherto unknown generalized effective medium approximation formulae are reported. Among them the analytical results (125) and (132) for the cylindrical cases are of greatest interest. The generalized theory recovers previously known anisotropic effective medium approximations in the literature as special cases. The practical applicability of the theoretical predictions has been tested by a simple comparison with experiment.
The primitive functions are [97] also for complex parameters d, e. For real parameters the same result can be seen more directly from the substitution s 2 = (u + (1/d 2 ))(u + 1)(u + e 2 ) in Eq. (85).