## Abstract

A material fluid element within a porous medium experiences deformations due to the disordered spatial distribution of the Darcy scale velocity field, caused by the heterogeneity of hydraulic conductivity. A physical consequence of this heterogeneity is the presence of localized kinematical features such as straining, shearing and vorticity in the fluid element. These kinematical features will influence the shape of solute clouds and their fate. Studies on the deformation of material surfaces highlighted the importance of stretching and shearing, whereas vorticity received so far less attention, though it determines folding, a deformation associated with the local rotation of the velocity field. We study vorticity in a three-dimensional porous formation exploring how its fluctuations are influenced by the spatial structure of the porous media, obtained by immersing spheroidal inclusions into a matrix of constant hydraulic conductivity. By comparing porous formations with the same spatial statistics, we analyse how vorticity is affected by the different shape and arrangement of inclusions, defined as the medium ‘microstructure’. We conclude that, as microstructure has a significant impact on vorticity fluctuations, the usual second-order statistical description of the conductivity field is unable to capture local deformations of the plume.

## 1. Introduction

Flows through porous media occur in a wide array of phenomena in natural and engineered systems. A distinctive property of a porous medium is the spatial variability of its hydraulic properties, which controls the rate at which chemical species react and disperse in geological formations [1–3], heat transfer [4,5] and water plant uptake [6,7]. The heterogeneous structure of the porous medium induces deformation of the fluid material elements which leads to spatially erratic flow paths. For these types of flow, the study of flow kinematics and topology is important given their strong role in controlling scalar dilution rates [8,9] and solute residence times [10,11]. For example, it has been shown that the occurrence of flow focusing (e.g. compression of fluid material elements) in porous formations enhances mass transfer [12–15].

This paper examines vorticity, the anti-symmetric component of the deformation tensor (see e.g. [16], ch. 2.3) in three-dimensional spatially heterogeneous porous media, at the Darcy scale. In porous media, vorticity has been studied less than the symmetric component of the deformation tensor, perhaps because it vanishes in a homogeneous isotropic hydraulic conductivity field, regardless flow dimensionality. On one hand, in a two-dimensional heterogeneous hydraulic conductivity field, vorticity is orthogonal to the plane containing the two components of the seepage velocity and, therefore, it cannot induce deformations of the material elements. On the other hand, in three-dimensional flow fields rotations around a point, which is the type of movement vorticity encapsulates, can potentially induce significant deformations if the hydraulic properties of the medium are spatially variable [17,18].

The statistical properties of vorticity in a heterogeneous random Darcy-scale velocity field have been investigated by Shvidler [19] and Kapoor [20]. However, their analyses are based on perturbation theory and, therefore, are limited to low heterogeneity of the hydraulic conductivity field. Within the context of numerical modelling of flow in geological formations, Mahani *et al.* [21] and Ashjari *et al.* [22] suggested to use vorticity as an indicator of heterogeneity and for improving coarse grid generation, as it accounts for both flow and permeability variations in a single measure. Furthermore, a detailed description of the vorticity field is fundamental for the computation of topological metrics, such as helicity [23] and the Okubo-Weiss parameter [24,25]. Recently, de Barros *et al.* [8] showed how enstrophy (i.e. vorticity squared) together with stretching and shearing deformation could be used to measure the mixing potential of the porous formation. Sposito [26], and more recently Chiogna *et al.* [27], studied vorticity in both locally isotropic and anisotropic conductivity fields. Several studies [28,29] showed how non-uniformity in the statistical properties of the hydraulic conductivity field can lead to spiralling flows, which affect solute mixing rates. Similar kinematics analysis was conducted by Raats [30].

Motivated by the aforementioned works, we study the vorticity field in a spatially heterogeneous three-dimensional porous formation. As opposed to previous contributions, we systematically investigate how the vorticity fluctuations, represented by the spatial mean of the vorticity's magnitude, are affected by medium heterogeneity. To achieve this objective, we make use of the multi-indicator model (MIM) of permeability structure [31,32], by which the porous formation is represented as an ensemble of inclusions of independent random hydraulic conductivity. This conceptualization is flexible and allows to obtain accurate solutions of the velocity field regardless the degree of heterogeneity. More precisely, the present approach allows consideration of any types of conductivity distribution, including broad probability density functions (pdfs) with infinite variance [33]. In particular, MIM provides a series expansion type solution, which reduces to an analytical formulation in the dilute case, i.e. when inclusions do not interact [32,34]. We emphasize that, generally in highly heterogeneous conductivity fields, solutions of the flow equations are very difficult to obtain, either analytically or numerically. For the sake of clarity, by means of this model it is possible to generate different types of media as a function of a few relevant parameters, which depend on the inclusion shape and its spatial arrangement, denoted in following as the medium ‘microstructure’. Such microstructure may for instance represent the orientation and spatial features of the layers that is often observed in aquifers. In this work, we are particularly interested in studying formations composed by different microstructures, but characterized by the same statistics of the spatial variability, such as the mean, the variance and the integral scales of the hydraulic conductivity field. The latter properties describe the medium ‘macrostructure’, i.e. a global, second-order-based description of the heterogeneous conductivity field. A qualitative example of such kind of media is shown in figure 1. However, this statistical description of the hydraulic conductivity field is truncated to the second-order, which may provide a lacking picture of the flow deformations, especially for high heterogeneity [35]. Thus, porous formations characterized by the same spatial statistics may show a very different flow behaviour. Motivated by these observations, our purpose is to fundamentally understand the impact of small-scale sedimentary structures on the vorticity field, as a measure of the rotational deformations of material elements caused by the hydraulic conductivity heterogeneity.

## 2. Mathematical framework

### (a) Conceptual model of the porous medium structure

With the goal of quantifying flow in heterogeneous porous media, we conceptualize the formation as a composite medium, made up of an ensemble of independent, non-overlapping *N* spheroidal inclusions of distinct hydraulic conductivity submerged into a homogeneous porous matrix of hydraulic conductivity *K*_{0}. With ** ζ**=(

*ζ*

_{1},

*ζ*

_{2},

*ζ*

_{3}) denoting the local Cartesian coordinate system with origin centred in the spheroid's centre of mass, the external surface of the spheroid is given by the following expressions

*a*and

*c*are the two main axes, respectively, the largest and the smallest, of the spheroid. The shape of the spheroid is characterized by the eccentricity

*φ*=

*c*/

*a*is the axis ratio. The inclusion volume varies with the axis ratio according to the following expressions:

*w*(

*φ*)=4

*πa*

^{3}

*φ*

^{2}/3 and 4

*πa*

^{3}

*φ*/3 for prolate and oblate spheroids, respectively [36], p. 131. For

*w*=4

*πa*

^{3}/3.

The inclusions are placed randomly within a domain of volume *W*, such that the volume fraction *n* occupied by the inclusions is equal to
*K*_{j} of the generic inclusion *j* (*j*=1,…,*N*) is either constant or randomly drawn from a lognormal pdf, depending on whether the medium is binary or unimodal, respectively. In the latter case, the hydraulic conductivities of the inclusions are independent and identically distributed random variables. However, we emphasize that the advantage of this method is that it allows the use of any conductivity distributions, regardless of the degree of heterogeneity [31,32]. The generic inclusion *j* is further characterized by the type of spheroid (oblate or prolate), the centroid location *e*_{j} and the orientation *l* and *m* identify the direction of the non-revolution axis with respect to the general reference system **x**=(*x*_{1},*x*_{2},*x*_{3}); the set of parameters is hereinafter denoted by the vector *θ*_{j}. A sketch of the porous media is shown in figure 2*a*, while the angles used to identify the orientation of the inclusions are illustrated in figure 2*b*. Based on these positions, the medium properties are fully described by the joint probability of the parameters *θ*_{j}. This conceptual model for the porous medium is similar to that employed in previous studies dealing with effective properties in composite media [37,38] and in flow and transport processes in complex subsurface environments [39,40].

The statistical properties of this composite medium are provided by Dagan *et al.* [41] for the simple case of spherical inclusions. The case of spheroidal inclusions is more complex, requiring a specific analysis of the correlation lengths as a function of the statistical properties of the parameters *θ*_{j}. However, an important payback of this higher complexity is in the flexibility the model offers in reproducing the correlation features of the hydraulic conductivity observed in the field [31], including statistical anisotropy [35], table 2.2. Such analysis is carried out in the following.

The hydraulic conductivity field *K*(**x**) of the spatially variable porous media can be expressed as follows:
*ξ*(**x**) is an indicator function which is equal to one when **x** is within the inclusion centred in *ξ* have been thoroughly discussed in [41,31]. Since *K*_{j} are independent and identically distributed random variables, the statistics of *K*(**x**) depends on the joint pdf of the parameters *θ*_{j}. Owing to the statistical independence of the *N* values of *K*_{j}, the ensemble mean of the hydraulic conductivity (2.4) assumes the following expression: 〈*K*(**x**)〉=〈*K*〉*n*+*K*_{0}(1−*n*). Similarly, the two-point covariance function *C*_{K}(**x**,**y**)=〈(*K*(**x**)−〈*K*〉)(*K*(**y**)−〈*K*〉)〉 is given by:
*χ*_{jk} is the overlap function representing the normalized intersection volume of two inclusions *j* and *k* shifted by **x**−**y**; the normalizing volume is the volume of the inclusions such that *χ*_{ij} varies between *χ*_{ij}=0 (no overlapping) and *χ*_{ij}=1 (complete superimposition) (see [41], section 2.1).

The above statistics have been derived under the assumption that the inclusions are far enough to neglect their mutual interference, i.e. the media is dilute (this point shall be retaken later), such that the terms of order *O*(*N*) have been retained in the calculation of (2.5) and the following condition applies: *j*≠*k*, which derives from the independency between the inclusions.

The integral scale along the *p* direction (with *p*=1,2,3) is defined as follows:
**r**_{p} the two-point separation vector along the direction with respect to which the integral scale is evaluated and *r*_{p}≡|**r**_{p}| is its module. According to expression (2.6), the directional integral scale depends on the mean of the overlap function *χ*_{jk}(**x**−**y**), which in turn depends on the distribution of shapes and orientation of the spheroids. The detailed procedure used for calculating the directional integral scales is provided in appendix A. The directional integral scales *I*_{p} (2.6) are generally employed as basic descriptors for the spatial correlation of the hydraulic conductivity in aquifers (see e.g. [35,40], ch. 4.6.4). However, in highly heterogeneous formations, the second-order description of the porous medium is not exhaustive and different behaviours can be observed in different formations with the same integral scales, as will be shown in the sequel with reference to vorticity.

### (b) Flow solution and vorticity

We consider steady-state flow of an incompressible fluid in an unbounded domain. The mean flow field is given by **U**=(*U*,0,0), where *U* denotes the mean seepage velocity aligned with the *x*_{1}-direction. The flow is divergence free and governed by the following mass conservation equation:
**V** is the Darcy-scale velocity, *ϕ* is the effective porosity (hereinafter assumed as constant and set, for simplicity, as unit) and *H* is the hydraulic head.

Solving equations (2.7) and (2.8), assuming an infinite domain and the spatially heterogeneous conductivity field (2.4), is a formidable task [31,34]. In order to grasp the main features of the vorticity field without resorting to computationally expensive numerical methods, we introduce the assumption that the media is dilute, i.e. *n*≪1. The dilute approach has a long tradition in fluid mechanics (see, e.g. [39,42,43]) and it is quite effective in deriving approximate analytical flow solutions in complex composite media, which are able to capture the main features of the velocity field and the associated vorticity field. Under such assumption, the inclusions are separated enough such that their interaction is weak and can be neglected.

This way, the velocity **V** can be written as the superposition of the velocity disturbance caused by each inclusion, as follows:
**v** is the fluctuation of the velocity field determined by an isolated inclusion *j* centred on *θ*_{j}, submerged in a porous matrix of conductivity *K*_{0}, subject to an imposed mean velocity **U** at infinity in a very large domain. Fully analytical solutions for **v** are available from the literature (e.g. [44], p. 427).

The vorticity ** Ω**(

**x**) of the velocity field (2.9) assumes the following expression:

*ω*_{j}is the vorticity field pertaining to the isolated inclusion

*j*. As (2.7) and (2.8) describe irrotational flows, the vorticity is confined at the surface of the inclusions, where the velocity typically displays a change in magnitude and direction, from the exterior to the interior of the inclusion.

Given that the mean velocity field is uniform and unidirectional, the mean vorticity in the system *Ω*(**x**)≡|** Ω**(

**x**)| indicates the magnitude of the vorticity field. After substituting (2.10) into (2.11), the SMVM assumes the following expression

*S*

_{j}is the spatial integral of

*ω*

_{j}≡|

*ω*_{j}| over

*W*(for further details, see appendix B). By taking the limit of (2.12) for

*I*

_{1}is the integral scale (2.6) in the longitudinal direction (

*p*=1).

Therefore, (2.14) is a fundamental relation, which we will use thoroughly in the present work, as it expresses a measure of the spatial variability of vorticity as a function of the statistical properties of the porous medium. The calculations leading to the expression (2.14) for spherical and spheroidal inclusions are detailed in appendix B.

## 3. Analysis

In this section, we investigate the influence of the medium structure on

### (a) Isotropic media

#### (i) Spherical inclusions

As a starting point, we consider first a binary isotropic porous medium, obtained by a random distribution of spherical inclusions with constant radius *A* and hydraulic conductivity *K* submerged in a matrix of uniform hydraulic conductivity *K*_{0}, such that the logconductivity contrast *y* is the same for all inclusions. The integral scale of the hydraulic conductivity for such medium is given by (4.4). The SMVM can be obtained analytically by substituting (4.17) into (2.14)

The behaviour of *y* is shown in figure 3 (green lines with circular markers). As expected, when the absolute value of *y* deviates from zero, namely when the hydraulic conductivity *K* of the inclusions deviates from *K*_{0}, *y*>0 than for *y*<0, as an effect of the larger deformation of the streamlines induced by flow focusing, a process that has been shown to develop only for *y*>0 [45,46].

Next, we consider a heterogeneous ensemble of inclusions in order to represent sedimentary formations of univariate lognormal conductivity. Following a customary procedure [31,32], the logarithm of conductivity contrast *y* is drawn from a normal pdf *f*_{y}, characterized by the mean

#### (ii) Spheroidal inclusions

As shown in appendix A, an isotropic medium can be obtained also with spheroidal inclusions if their two non-revolution axes are oriented randomly according to the uniform pdf (4.2), where both *α*_{l} and *α*_{m} are bounded between 0 and *π*/2. For spheroidal inclusions with uniform logconductivity contrast *y* (bimodal media) and a given axis ratio *φ*, the SMVM assumes the following expression
*I*_{1} and volume fraction *n*; hence the media made up of either spheres or spheroids mainly differ in their microstructure, i.e. the particular arrangement of spheroids and their spatial orientation.

Figure 3*a*,*b* show *y*<0 and *y*>0, respectively, and for both prolate and oblate inclusions with the following axis ratios: *φ*=1.00, 0.10 and 0.01. For both high and low contrasts, *φ* decreases, changing the inclusion shape from spherical to flattened oblate or elongated prolate spheroids. Considering the same value of *φ*, for low conductivity inclusions (*y*<0) (figure 3*a*), oblate spheroids always lead to a SMVM larger than for prolate spheroids. For high conductive inclusions (*y*>0) (figure 3*b*), a similar behaviour is observed for small to moderate contrasts, but a crossover occurs at values of *y* that increase as *φ* reduces, such that prolate inclusions yield the larger SMVM for large contrasts. In addition, the disparity between the upper bounds of SMMV increases as *φ* reduces. Indeed, the ratio between the asymptotic values of *φ*=0.1 and 0.01, respectively, and to 30 and 1800 for prolate inclusions with the same *φ* values. As can be seen, the axis ratio *φ*, which is the principal parameter describing the medium microstructure, plays a major role in controlling *φ* yields to: (i) an increase or a reduction of the single inclusion's SMVM, *S*(*α*_{l},*α*_{m};*κ*,*φ*), depending on whether the inclusions are highly (*y*>0) or low conductive (*y*<0), respectively. Regarding to the spheroid type, the higher *S*(*α*_{l},*α*_{m};*κ*,*φ*) is yielded by oblate inclusions when *y*<0 and by prolate inclusions when *y*>0; (ii) a reduction of the longitudinal integral scale *I*_{1}(*π*/2,*π*/2;*φ*) and (iii) a reduction of the inclusion volume *w*(*φ*). In particular, a reduction of *w*(*φ*) is balanced by an increase of the number of inclusions, within a given volume fraction (2.3). For instance, we observe in figure 1 that increasing the number of inclusions contained within a given volume enhances the disorder in the distribution of the hydraulic conductivity and, consequently, enlarges the deformation of the streamlines. These results are consistent with the significant deformations of the material surfaces already observed for highly conductive prolate inclusions by Janković *et al.* [47].

The case of a unimodal heterogeneous porous media with normally distributed *y* is shown in figure 4. For a relatively large axis ratio, i.e. *φ*=0.1, and weakly to moderate heterogeneity, the largest *φ*=0.01) show a similar behaviour, though with larger relative differences, at small to intermediate

### (b) Anisotropic media

Anisotropic conductivity fields can be obtained by limiting the range of random rotation of the spheroid's non-revolution axes to angles smaller than *π*/2. Rotations are chosen in such a way to obtain the desired statistical axisymmetric anisotropy of the hydraulic conductivity, *f*=*I*_{3}/*I* with *I*=*I*_{1}=*I*_{2} and *I*_{3}<*I*, according to (4.12). Details of the calculations are provided in the appendix A. Figure 5 shows *φ*=0.001 and the following anisotropy ratios: *f*=1.00, 0.10, 0.01 and 0.001 for oblate spheroids (lines with warm colours) and *f*=1.00, 0.50 and 0.25 for prolate spheroids (lines with cold colours).

Notice that for oblate inclusions, the lowest anisotropy ratio *f*_{min}=*I*_{3}/*I*_{1}=*φ*=0.001 is reached when all the inclusions are with the longest axis aligned to the main flow direction; on the other hand, for prolate inclusions the smallest anisotropy ratio, *x*_{3} axis, according to equation (4.13). We found that the statistical anisotropy impacts significantly *f* reduces; these results are in agreement with previous investigations [18,28,29].

Similarly to the isotropic case, the microstructure has a different impact on SMVM, depending on the conductivity contrast: at low to moderate contrasts the largest *φ*), while prolate inclusions overwhelm oblate inclusions at large contrasts. The crossover of prolate with respect to oblate inclusions at relatively high heterogeneity is due to the intense flow focusing induced by highly conductive prolate inclusions. Figure 6 shows *φ*=0.001 and the following anisotropy ratios: *f*=1.00, 0.10, 0.01 and 0.0001 for oblate spheroids (lines with warm colours) and *f*=1.00, 0.50 and 0.25 for prolate spheroids (lines with cold colours). SMVM grows with *f*. As for the binary case,

## 4. Summary

In this work, we analysed the role of the hydraulic conductivity structure of heterogeneous porous formations in developing vorticity, which is the controlling factor of folding, a deformation of the material surface that, together with stretching and shearing, determines mixing in natural aquifers. More specifically, our interest was to explore how the spatial mean of the vorticity's magnitude, *y*, i.e. the logarithm of the contrast between the hydraulic conductivity of the inclusions and the underlying matrix, leading to a binary or unimodal distribution, respectively. Then different microstructures are obtained by using both oblate or prolate spheroids and changing *φ*, the ratio between the two main axes. The vorticity's magnitude is normalized with respect to the ratio between the mean velocity and the longitudinal logconductivity integral scale, such as to compare cases with the same ‘macrostructure’, i.e. the same large-scale variability as described by second-order statistics. Large-scale structure of the hydraulic conductivity is therefore controlled by the distribution of *y*, while microstructure is dictated by the shape of the inclusions and their orientation with respect to the mean flow.

A first analysis, conducted with binary media, showed that *y*, approaching asymptotic values which depend on the microstructure. The more evident effect of microstructure is for highly conductive inclusions (i.e. *y*>0), which converge to the asymptotic values slower than the low-conductive inclusions (i.e. *y*<0). This can be explained with the larger deformations of the streamlines caused by the focusing effect produced by high-conductivity inclusions, an effect that is not observed when *y*<0. Additional analysis indicate that diminishing *φ* leads to a vorticity field with higher *n*, the reduction of the inclusion volume, caused by the reduction of *φ*, is balanced by a larger number of the inclusions. As a consequence, streamlines are more contorted, thereby leading to a larger *f* reduces and the field becomes more anisotropic. Notice that *et al.* [47], fig. 8.

Summarizing, we demonstrate that different microstructures, which may represent the types of layering occurring in natural porous formations, may significantly affect the vorticity distribution and the kinematical features of the flow field. Hence we emphasize that the second-order approximation, which is commonly adopted in studies dealing with random porous media, may be not adequate to picture the local deformations of the plume, weakening our ability to assess dilution and reaction rates of pollutants, particularly in highly heterogeneous porous media.

## Authors' contributions

All authors conceived the study and wrote the paper, A.F. and M.D.D. developed the analytical expressions of vorticity and M.D.D. performed the simulations.

## Competing interests

We declare we have no competing interests linked to this study.

## Funding

This research was partially supported by the Italian Ministry of Public Instruction, University and Research through the project PRIN 2010-2011 (Innovative methods for water resources management under hydro-climatic uncertainty scenarios, prot. 2010JHF437).

## Acknowledgements

We are grateful to G. Dagan for stimulating discussions on the subject matter. Furthermore, the first author gratefully acknowledges the University of Trento for the financial support provided by the PhD scholarship, C. B. Rizzo for the editorial assistance and the kind hospitality of the Sonny Astani Dept. of Civil and Environmental Engineering at the University of Southern California during her stay in the USA.

## Appendix A. Directional integral scales

By substituting expression (2.5) into (2.6), with the two-point separation vectors oriented along the axes of the coordinate system, we obtain the following expressions for the three main integral scales
*χ* is the overlapping volume of two inclusions, with the same shape and orientation, translated by −*r* in the same direction along which the directional integral scale is evaluated. *χ* is normalized with respect to the volume of the inclusions such that it assumes values bounded between 0 and 1 [40]. The two-point separation distance *r* along the main axes directions depends on the angles ** α**=(

*α*

_{l},

*α*

_{m}), identifying the orientation of the revolution axis of the spheroid with respect to the reference system (

*l*=2,

*m*=3 for prolate and

*l*=1,

*m*=2 for oblate spheroids, respectively), and

*p*is the probability density function of the rotation angles. Hereafter we will consider a uniform pdf of the rotation angles in the range [−

*γ*

_{l},

*γ*

_{l}] and [−

*γ*

_{m},

*γ*

_{m}] such that:

**(a) Spherical inclusions**

For spherical inclusions of radius *A*, the overlapping normalized volume is independent on the direction and assumes the following expression [40]:

**(b) Spheroidal inclusions**

In order to calculate the integral scales for spheroidal inclusions, the local reference system ** ζ**=(

*ζ*

_{1},

*ζ*

_{2},

*ζ*

_{3}), centered in the spheroid's centre of mass, is deformed by (i) a rotation by the angles

**=(**

*α**α*

_{l},

*α*

_{m}), so that the spheroid's major axes is aligned with the

*ζ*

_{1}coordinate, (ii) a magnification of the minor axis, so that the transformed spheroid becomes equal to a sphere of radius

*A*. Considering the above transformations, the integral scales for the spheroid are given by (4.4):

*A*

_{i}are the spheroid's axis, depending on the axis ratio

*φ*and the spheroid inclination

**, in the transformed reference system and**

*α**d*

_{i}are determined imposing the equivalence

*A*

_{i}=

*A*, and assume the following expressions:

The expression (4.5) is for a single spheroidal (prolate or oblate) inclusion inclined by the angles ** α**=(

*α*

_{l},

*α*

_{m}) with respect to the global reference system

**x**=(

*x*

_{1},

*x*

_{2},

*x*

_{3}). For a large number of inclusions with varying inclination, within the hypothesis of a dilute systems, the integral scales can be computed as the ensemble average of expressions (4.5) with the pdf of the inclination angles given by (4.2).

In the case of statistical axysimmetric anisotropy, with anisotropy ratio *f* as defined in §3*b*, the integration limits *γ*_{l} and *γ*_{m} in (4.1) are selected by solving the following system of two equations
*γ*_{l} and *γ*_{m} obtained from the solution of (4.12) are then substituted into equations (4.1) to obtain the directional integral scales *I*_{i}. According to this model, for oblate spheroids, the highest anisotropy is *f*=*I*_{3}/*I*_{1}=*φ*, which is obtained for *γ*_{l}=*γ*_{m}=0 and all the inclusions are with the longest axis aligned to the main flow direction. On the other hand, for prolate spheroids, the largest anisotropy is obtained by rotating randomly the inclusions around the *x*_{3}-axis between *α*_{3}=−*π*/2 and *π*/2; the resulting anisotropy ratio is

## Appendix B. Calculation of S ¯ ( θ ) for an isolated inclusion

**V**^{ex} and **V**^{in}, are calculated on the inclusion's boundary ∂*w*. Normal components are conserved, whereas tangential components are discontinuous across the interface. (ii) We guess that the transition between exterior and interior velocity occurs linearly along a *d* thick shell surrounding the inclusion. (iii) Into the shell, velocity field is described by a local cartesian system (*x*,*y*,*z*) centred in a point on the inclusion's boundary, with *z* normal to the surface: **V**(*x*,*y*,*z*)=**V**^{in}|_{∂w}−(**V**^{in}|_{∂w}−*V*^{ex}|_{∂w})*z*/*d*. (iv) As local velocity is a function of *z* solely, the curl of velocity into the shell is constant and depends only on the velocity gradient that occurs along the shell (*V*^{in}|_{∂w}−*V*^{ex}|_{∂w})/*d*. (v) Finally, the curl's magnitude is integrated on the shell volume: first normally to the surface of the shell and then along it. The former integration removes *d*, thereby eliminating the need to compute the limit for

Following all this step, the spatial mean of the vorticity's magnitude for a single inclusion becomes:

**(a) Spherical inclusion**

For a single spherical inclusion centred at the origin in a dilute system (*n*≪1), the solution for the potential *ϕ* is taken from [41]
*A* is the sphere radius and *β*=(1−*e*^{y})/(1+*e*^{y}/2). The velocity field is obtained computing the gradient of (4.15), using spherical coordinates system *r*=*A*), interior and exterior velocities are

Therefore, the spatial mean of vorticity's magnitude for a single spherical inclusion depends on the mean velocity *U*, on the radius *A* and on the value of logconductivity contrast *y* (which is represented by *β*).

**(b) Spheroidal inclusion**

For a single prolate or oblate spheroidal inclusion, the potential *ϕ* assumes the expression provided by Carslaw & Jaegger [44], p. 427, in the study of the conduction of heat in a spheroid of thermal conductivity *K* plunged in to a medium of constant thermal conductivity *K*_{0}. These expressions, based on the Maxwell's solutions, are omitted here for the sake of conciseness. Unlike the case of a spherical inclusion, the integral (4.14) does not admit analytical treatment, such that it has been computed by means of a suitable numerical quadrature.

- Received October 23, 2015.
- Accepted February 19, 2016.

- © 2016 The Author(s)

Published by the Royal Society. All rights reserved.