Analytical solutions for deformable elliptical inclusions in general shear

SUMMARY Using Muskhelishvili’s method, we present closed-form analytical solutions for an isolated elliptical inclusion in general shear far-field flows. The inclusion is either perfectly bonded to the matrix or, as in the case of a circular inclusion, to a possible intermediate layer. The solutions are valid for incompressible all-elastic or all-viscous systems. The actual values of the shear modulus or viscosity in the inclusion, mantle and matrix can be different and no limits are imposed on the possible material property contrasts. The solutions presented are complete 2-D solutions and the parameters that can be analysed include all kinematic (e.g. stream functions, velocities, strain rates, strains) and dynamic parameters (e.g. pressure, maximum shear stress, etc.). We refrain from giving the tedious derivation of the presented solutions, and instead we focus on how to use the solutions, how to extract the parameters of interest and how to apply and verify them. In order to demonstrate the usefulness of Muskhelishvili’s method for slow viscous flow problems, we apply our results to a clast in a shear zone and obtain important new insights, even for simple, circular inclusions. Another important application is the benchmarking of numerical codes for which the presented solutions are most suitable due to the infinite range of viscosity contrasts and the strong local gradients of properties and results. To stimulate a broader use of Muskhelishvili’s method, all solutions are implemented in MATLAB and downloadable from the web.

. Setup of circular inclusion problem. A clast with radius r c is embedded in a matrix and subjected to combined pure and simple shear far-field flows.
In the case of a mantled circular inclusion, a layer of constant thickness is introduced between clast and matrix. The layer-matrix interface is defined through radius r l . The origin of all coordinate systems is chosen to be the centre of the clast. The position of a point P in the z-plane can be expressed by three different coordinate systems, P(z p ), P(x p , y p ), P(r p , θ p ). Note that the sketched box size is not to scale; the presented analytical solutions are based on the assumption that the box boundaries are far from the clast. C 2003 RAS, GJI, 155, 269-288 Downloaded from https://academic.oup.com/gji/article-abstract/155/1/269/713923 by guest on 28 July 2018 geological system that can hamper the pressure-temperature path interpretation and hence, the conversion of pressures into burial depth (e.g. Tenczer et al. 2001). Another shortcoming of the mentioned solutions is that they are strictly only valid for infinite viscosity contrasts between the clast and matrix. Since finite viscosity contrasts and slipping interfaces are more likely to be the norm than the exception, a corresponding analytical solution for the complete kinematics and dynamics is necessary. The similarity between the mantled clast and the slipping interface clast is that at the limit of vanishing mantle viscosity and small mantle thickness, the normal traction and velocity are continuous through the interface, but the shear tractions vanish and the tangential velocities are discontinuous, characteristic of a slipping interface. Thus, the slipping interface clast is an end-member case of the mantled clast.

Numerical code benchmarking
The advantages of the presented analytical solutions is that they are complete 2-D closed-form solutions and come with a set of rules on how to extract any desired kinematic or dynamic parameter. Yet, cases for which analytical solutions can be found are rare and generally asymptotic methods (e.g. Barenblatt 1996) combined with numerical models must be used in addition or as an alternative. Nevertheless, analytical solutions retain their importance because numerical codes must be tested thoroughly (see Pelletier et al. 1989, for possible numerical problems with simple incompressible viscous flows). The majority of analytical benchmark tests in use are restricted to essentially 1-D viscosity profiles (Moresi et al. 1996). 'Essentially 1-D' refers to the category of tests that are based on linear stability analysis such as 2-D folding (e.g. Biot 1961;Fletcher 1974) or diapirism (Chandrasekhar 1961). Mathematically, the viscosity contrast does not affect the left-hand side in these problems, i.e. the stiffness matrix of the system, but only contributes to the right-hand side of the system of linear equations generated by the discretization process. These problems do not account for any interaction between different harmonics (in frequency space) and are therefore only partly suitable for 2-D benchmarking.
The solutions presented here implement genuine 2-D problems. The complete solution is described using simple polynomials for complex potentials. The entire viscosity contrast range is covered and strong local gradients of properties and solutions may be present (such as within the rim), features that are crucial for proper code benchmarking.

Solution derivation, implementation and availability
The interest in isolated lubricated or perfectly bonded particles is not restricted to the geological community but is, in fact, inherent to many other fields of science. The importance of the behaviour of composites has stimulated researchers in recent years to derive the analytical solution for problems similar to the one posed here. As expected, the methods used are those of Eshelby (1957Eshelby ( , 1959 and Muskhelishvili. Relevant references include Mura (1987), Stagni (1991), Furuhashi et al. (1992), Huang et al. (1993), Gao (1995), Ru & Schiavone (1997), Shen et al. (2001), and the basic building blocks of the solutions can already be found in the books of Muskhelishvili (1953) and Savin (1961).
Common findings of the more recent works are that: (1) the Eshelby conjecture (e.g. Mura 2000) does not hold for imperfect bonding between matrix and clast and (2) the lubricated ellipse does not have a closed-form solution. The first finding simply means that the stress and strain (rate) components inside the inclusion cannot be described by a single value, which is the case for perfect bonding between clast and matrix. The second finding means that the solution is an infinitely long series (Ru 1999;Shen et al. 2001). For this reason we refrain from presenting the mantled elliptical inclusion here. Instead we concentrate on the cases for which finite series solutions can be found.
Our solutions fulfil mass and force balance, rheological equations and far-field boundary conditions. They are obtained by matching the tractions and the velocities through the interfaces and solving the linear system of equations for the coefficients of the polynomial representations of complex potentials. In the same manner it is possible to verify the presented solutions.
All solutions given have been implemented in MATLAB and can be downloaded from the web (http://e-collection.ethbib.ethz.ch/show? type = bericht & nr = 188) or from the authors upon request. The figures in this paper are reproducible in colour with these MATLAB scripts.

B R I E F R E V I E W O F M U S K H E L I S H V I L I ' S M E T H O D
An in-depth introduction to Muskhelishvili's method is beyond the scope of this paper. We simply provide the most important equations and describe how the various parameters of interest can be extracted from the complex potential solutions. Interested readers may consult Muskhelishvili's original work, Lu (1995), or Jaeger & Cook (1979 for further documentation.

Basic set of equations
Muskhelishvili's method makes use of the fact that problems in two dimensions can be conveniently expressed in terms of a complex coordinate z. As shown in Fig. 1 the coordinates of a point in the complex plane can be expressed in different coordinate notations: where i = √ −1, x and y are the usual Cartesian coordinates, and r and θ are the radius and the orientation angle of the polar coordinate system, respectively. C 2003 RAS, GJI, 155, 269-288 The basis of Muskhelishvili's method is that the bi-harmonic equation, which describes the 2-D plane stress or plane strain elasticity problem, has a general solution that can be expressed in terms of two complex functions, φ(z) and ψ(z). The conditions imposed on φ(z) and ψ(z) are that they must conform with the applied boundary conditions and be analytic, i.e. satisfy the Cauchy-Riemann equations (e.g. Jaeger & Cook 1979). In simple terms, this means that complex potentials must be 'normal' functions of z, such as sin (z), ln (z), z 2 , but notz, (z) or (z). The overbar denotes conjugation, signifies the real part and the imaginary part.
Since we advocate a more widespread use of the method to slow flow problems, we provide here the three fundamental equations of the Muskhelishvili method for slow, incompressible, viscous flow in plane strain: where σ xx , σ yy and σ xy are the components of the total stress tensor, v x and v y are the horizontal and vertical velocities, µ is the viscosity of the material for which φ(z) and ψ(z) are valid, prime and double prime denote the first and second derivatives versus z. Once the analytical expressions of φ(z) and ψ(z) are obtained, stresses, velocities and a variety of other parameters can be evaluated. For example, the pressure perturbation (p) is obtained through eq. (2) as The force balance equations determine the pressure only up to a constant. Therefore, an arbitrary (lithostatic) value may be added without any influence on the result. The sign convention used is that compressive pressure perturbations are positive. Another useful parameter is the effective or maximum shear stress (τ ), which is calculated as (e.g. Ranalli 1995) The Cartesian components of the velocity, v x and v y , are obtained from eq. (4) by taking the real and imaginary parts, respectively. Another practical parameter for the analysis of 2-D flows is the stream function, (Turcotte & Schubert 1982). The contour lines of are used to visualize the flow of individual particles in the steady state. The Cartesian definition of is Hence, can be obtained by integrating either v x or v y .

Far-field flow expressions in complex potentials
The kinematic boundary conditions used in this work are the pure shear (ps) strain rate,ε, and/or the simple shear (ss) strain rate. In terms of complex potentials, we write: which through eq. (4) produce v x + iv y = zε using the identity given in eq. (1), the pure shear flow field is obtained, Simple shear

General boundary conditions
The addition of these shear flows yields the general, combined pure and simple shear (gps) boundary conditions: The distinction of simple and pure shear flow is somewhat redundant as at any instant one can be translated into the other through a coordinate system rotation. However, the provision of both facilitates the analysis of general shear flows.

C I RC U L A R I N C L U S I O N Solution
The circular inclusion in a matrix of different viscosity subject to the general boundary conditions can be solved with Muskhelishvili's method and a detailed account of this problem has been given in Jaeger & Cook (1979). Naturally, we need two sets of φ(z) and ψ(z) to describe the result. One set, φ(z) c and ψ(z) c , describes the result within the inclusion/clast and is valid from the origin, chosen to be the centre of the clast, to the clast radius, r c . The second set of complex potentials, φ(z) m and ψ(z) m , determines the solution in the matrix. Subscripts 'c' and 'm' are used to distinguish clast and matrix values, respectively. Clast

Clast kinematics
With this solution in terms of a complex potential solution, it is possible to derive useful expressions such as the complete kinematics of the clast by plugging φ(z) c and ψ(z) c into the velocity expression (eq. 4): This expression is valid for all possible viscosity contrasts between clast and matrix and for arbitrary combinations of simple and pure shear. It allows us to study the rotational behaviour of a circular clast in a shear zone. Jeffery (1922) has shown that the rigid inclusion rotates with a rate that is half of the applied shear rate. In order to reproduce this result, we need to derive the rotation rate,ω, that is inherent to eq. (26). In Cartesian coordinates, the rotation rate is defined as and is closely related to the shear strain rate,ε xy , that iṡ Hence, in simple shear only (ε = 0) the clast rotation rate iṡ which is identical to Jeffery's result. If we apply a top to the right (positive) shear strain rate, then the clast rotates with the simple shear flow in a clockwise sense, at a rate that is half the applied shear rate. However, the rotation rate that we derived is not only valid for rigid clasts, but for arbitrary viscosity contrasts between clast and matrix. Therefore, the factor of 2 difference between the applied shear rate and the rotation rate persists independently of the actual viscosities. However, only the kinematics of the infinitely rigid inclusion are reducible to a rigid-body rotation. All other circular clasts will show some shear deformation that iṡ As predictedε xy → 0 for µ c → ∞ and the maximum shear rate is obtained for the infinitely weak inclusion, µ c → 0. When µ c = µ m , the shear strain rate in the clast is equal to the far-field value in the matrix. It is noteworthy that a clast that is only 100 times more competent than the matrix will already exhibit a shear deformation that is smaller than 2 per cent of the matrix value.

Clast dynamics
One of the driving forces of phase transitions and metamorphic reactions is pressure. The pressure inside the clast is obtained by plugging φ(z) c into eq. (5), which yields Eq. (31) holds irrespective of the boundary conditions and viscosities-the clast itself remains at the background (lithostatic) pressure. The second important driving force is the differential or maximum shear stress, τ c . Analogous to the pressure, τ c inside the inclusion has the property that it can be described with a single constant value In the case of an infinitely weak inclusion, the maximum shear stress vanishes. In contrast, for the rigid inclusion, τ c is a function of the kinematic boundary conditions and the viscosity of the matrix This value is likely to be characteristic for most competent clasts since a viscosity contrast between clast and matrix of 10:1, already yields a maximum shear stress value that deviates less than 10 per cent from the µ c → ∞ case.

Matrix dynamics
The matrix values of p and τ do not show the property of single, constant values, but are by no means more complicated to derive. Using the given complex potentials of the matrix, the expression for pressure in polar coordinates is If we are interested in the matrix pressure on the clast-matrix boundary we simply set r = r c . We can immediately see that in simple as well as pure shear there are two pressure maxima and minima along an entire clast circumference. The amplitudes of these extrema depend on the viscosities and the applied far-field flow, but not on the clast size. This is also the case for all other stress and strain rate components and contradicts the proposition of Passchier & Simpson (1986) that the amplitude of the driving forces of clast recrystallization should decrease with decreasing clast size. It is also important to realize that pressure, in contrast to for example traction components, is not required to be continuous through an interface ( p c = 0 within the inclusion). As it was the case for the differential stresses in the clast (τ c ), eq. (34) shows that the infinite viscosity contrast pressure limit is approached rapidly and therefore this limit may be taken to be representative for the majority of natural clasts. The value of this infinite-pressure limit (in simple shear) is obtained by setting µ c → ∞ and θ = 3/4π , which results in Therefore, the maximum overpressure in the matrix is twice the applied far-field stress (µ mγ ). Eq. (34) is equally applicable to competent as well as to weak inclusions (µ c < µ m ). Thus we are able to plot the matrix pressures extrema at the interface and θ = 3/4π for the entire range of µ c /µ m (Fig. 2). The behaviour of the weak inclusion is exactly the inverse of the competent inclusion, i.e. the loci of compression become extensive and vice versa (cf. Figs 10a and b). This is interesting, because the applied simple shear far-field flow is still the same, i.e. the top left and the bottom right quarter (relative to the clast) are still 'streamed' at by the background component of the flow. However, if the clast is weaker than the matrix these quadrants go into relative extension, while the originally extensive quadrants become compressional (cf. fig. 6b in Tenczer et al. 2001). This is characteristic for most dynamic and kinematic parameters. The relatively simple eq. (34) contains the complete pressure information of the entire matrix as displayed for simple shear in Fig. 10(a). The size and the shape of the overpressure and underpressure regions can be captured and hence the estimated dimensions of possible pressure shadows. The pressure field, as expected, shows perfect point symmetry around the clast centre. Using a different method Masuda & Ando (1988) could not obtain this perfect point symmetry although their analytical series solution had 24 terms.

C I RC U L A R I N C L U S I O N W I T H A R I M Solution
The circular inclusion with a rim is a three-phase system where a zone of constant thickness (subscript 'l' for 'layer' in the following) and possibly different viscosity is introduced between the clast and the matrix. Accordingly, the solution consists now of three different sets of φ(z) and ψ(z) for clast, matrix and layer. Note that the solutions in clast and matrix are also altered due to the presence of the rim.
The additional parameters introduced are the viscosity of the layer, µ l , and the radius of the layer-matrix interface (Fig. 1). For geometrical reasons the condition imposed on r l is r l ≥ r c . If r l = r c the rim has zero thickness and the previous solution is recovered. In order to reduce the number of parameters involved we normalize by µ m and r c . Therefore, the viscosities present areμ c andμ l and the remaining radius is r l , where the tilde denotes normalization.
As a result of the increasing complexity of the geometry and number of parameters, the solution for the mantled circular inclusion is significantly more involved in terms of coefficients, but not in terms of the form of the solution. Clast Layer The coefficients Q 1 to Q 8 depend only onμ c ,μ l andr l , are therefore real, and are given in the Appendix. The general form of the solution in terms of which powers of z are present can be explained as follows. Inside the clast only positive powers of z are admissible because otherwise the functions would not be analytic since division by zero would occur in the centre of the clast (z = 0). On the other hand, in the matrix all values of stress, strain rate and velocities, and hence φ(z) and ψ(z), must decay towards the background/far-field state with increasing distance from the layer/clast. Therefore, the solution in the matrix consists of the far-field flow and negative powers of z. The complex potentials within the layer have the task of matching the matrix and the clast and hence the polynomials contain both, negative and positive powers of z.

Clast kinematics
We have seen that the circular clast with perfect, direct bonding to the matrix rotates in simple shear with a rate that is half the applied shear rate, irrespective of µ c /µ m . The question that arises is what influence an intermediate layer, strong or weak, between clast and matrix has on the rotation rate. Evaluating the Cartesian velocity components within the clast we obtain v x r c = + 1 2ỹ We can see that the first term in the simple shear part is a rigid-body rotation at a rate that is half the applied shear rate. Settingε = 0 and evaluating the expression for the rotation rate (eq. 27) results iṅ which confirms the above observation. However, clearly the Eshelby conjecture is violated since the rotation rate is not constant within the clast. Yet, for the infinitely rigid clast the second term of eq. (44) vanishes. This result is noteworthy since it implies that the rotation of a circular, rigid clast is independent of a strong (µ l > µ m ) or very weak/lubricant rim (µ l µ m ). We conclude here our investigation concerning the rotation rates of circular clasts: the rigid circular clast always rotates withω = −γ /2, irrespective of the matrix viscosity and the presence of additional interfaces/layers between clast and matrix. If the clast is not rigid the presence of a layer between the clast and the matrix disturbs this equality. Yet, if the clast is directly and perfectly coupled to the matrix it will always rotate withω = −γ /2, irrespective of the viscosity contrast between clast and matrix.

Dynamics
The exceptional characteristics of the circular clast embedded in a matrix and subject to uniform far-field boundary conditions are that the inside of the clast shows constant dynamic values and, in particular, no pressure perturbation is generated in simple shear. Plugging the φ(z) of the three different phases into the pressure equation the following expressions are obtained for simple shear only: Eq. (45) generally indicates that the property of zero-pressure perturbation in the clast is lost. It is interesting to observe that the presence of the rim has the effect of synchronizing the clast and the matrix. Both eqs (45) and (47) exhibit two minima and two maxima around the clast circumference at exactly the same locations. In the matrix, this was already obtained for the case of perfect direct bonding between circular clast and matrix. Indeed, eqs (47) and (34) are very similar and naturally are identical if the viscosity of the layer is equal to the matrix value.
Although the property of zero-pressure perturbation within the clast is lost, certain limits still exhibit this characteristic.
(1) When the viscosity of the layer is equal to the matrix value (μ l = 1) then Q 1 is zero and consequently p c = 0. (2) If the layer becomes lubricant (μ l → 0) Q 1 also vanishes and the clast pressure is lithostatic.
Plotting the complete 2-D pressure perturbation, stream function and maximum shear stress fields is straightforward. The only information required areγ ,ε,μ c ,μ l andr l . An example is given in Fig. 3

Conformal transformations
The genuine strength of Muskhelishvili's method is that conformal transformations can be applied. Conformal transformation means a coordinate transformation that is analytic, i.e. fulfils the Cauchy-Riemann equations. The geometrical meaning of conformal transformation is that the mapping is unique and that both amplitude and sense of angles are preserved. The reason why conformal transformations are useful is because they allow solution finding for complex physical problems (z-plane) in geometrically simpler image planes (ζ ). The most famous usage of conformal solution transformation may be by due to N.E. Joukowski who, in 1906, calculated the lift of an airfoil. He reduced the problem to calculating the flow around a rigid cylinder in the image plane and then, using what is now called 'Joukowski transform', translated his solution to the physical domain, where the cylinder becomes an airfoil (this class of airfoils now being called 'Joukowski airfoils').
The Joukowski transform is not only suitable for studying airplane wings, but is also appropriate for the study of elliptical inclusion. The most general definition of the Joukowski transform is However, for our study R and m can both be set to 1, since for elliptical shapes they are not needed. The characteristic effect of the Joukowski transform is explained in Fig. 4. Certain off-centre circles in the ζ -plane become airfoil-like shapes in z. Since we are interested in elliptical inclusions in z, we restrict the possible circles in ζ to the class for which the centre coincides with the coordinate system origin. As a result of the condition of uniqueness of the mapping we have to restrict the possible circle radii, ρ, to ρ ≥ 1. The circle with ρ = 1 becomes a slit in z and has an aspect ratio, t = a/b, which is infinite. The relationship between ρ and t is inverse, i.e. the aspect ratio in z decreases as ρ increases in ζ : Hence, for a given physical problem with an elliptical inclusion of aspect ratio t we can set up the corresponding problem in z, where we have the matrix from infinity to ρ c and the inclusion from ρ c to 1. Therefore, the inclusion is a ring in ζ (Fig. 4).

Far-field flow conditions
The far-field flow conditions in the physical domain are still pure shear and/or simple shear. Since we solve the problem in the image plane ζ , we must translate the boundary conditions from z to ζ . This is done by replacing z with the expression of the Joukowski transform. Given that we are dealing with a non-circular inclusion the angle (α) between the boundary condition flow and the long axis of the ellipse must also be considered (Fig. 5). The definition of α used is that it is measured from the horizontal of the boundary conditions and positive values are counter-clockwise. For this system we can write the combined pure and simple shear boundary conditions with the inclination α as Setting α = 0 restores the general horizontal boundary conditions (eqs 19 and 20).

Solution
The solution for the elliptical inclusion with inclined arbitrary combinations of pure and simple shear is where B contains rotated boundary condition terms and B 1 to B 5 are real-valued combinations ofμ c and ρ, As before we have used µ m to normalize the viscosities, and the inside radius of the inclusion ring in ζ to normalize all length parameters, hence the tildes.

Basic set of equations
The solution presented above is written in terms of the image plane coordinates, ζ . The original basic set of equations at the beginning of this paper was given for the physical plane z. To allow the analysis of the solutions we must translate the basic set of equations into the ζ -plane, taking the specific conformal mapping into account. Writing φ and ψ instead of φ(ζ (z)) and ψ(ζ (z)) the basic set of equations under Joukowski mapping becomes: .

Clast kinematics
Based on Muskhelishvili's method, Ghosh & Ramberg (1976) have derived the kinematic behaviour of the rigid elliptical clast in combined pure and simple shear. With the solution provided here it is possible to obtain the expression for the kinematic behaviour of any kind of elliptical inclusion, whether competent or weak, in combined pure and simple shear. Analysing the complex potential expressions of the clast it is evident that it is possible to rewrite them directly for the z-domain. Both φ(ζ (z)) c and ψ(ζ (z)) c are functions of (ζ + 1/ζ ), which is the definition of the used Joukowski transform (eq. 48). Therefore, we can replace (ζ + 1/ζ ) by z and obtain expressions that are valid directly in z. This greatly simplifies the analysis, since coordinate transformations are no longer employed and consequently the simpler z-set of the basic equations can be used. The resulting expression for the rigid body rotation rate iṡ This expression successfully reproduces Ghosh & Ramberg's (μ c → ∞) case and can be used to explore the entire field ofμ c − t. We can, for example, setε = 0 and varyμ c and t as is done in Fig. 6 for different aspect ratios, t. It can be observed that if the clast aspect ratio is big enough and the viscosity is less than the matrix value, a field of back rotation comes into existence (Fig. 7b), which is limited by the intersections with theω = 0 line. If the inclusion is inclined with an angle corresponding to this field, it will rotate against the applied shear sense towards the lower intersection point. Since this field of back rotation only exists for weak inclusions, it must be realized that the corresponding shear strain rates are not negligible and may overprint the back rotation. In the case of the infinitely weak inclusion, the minimum required aspect ratio for which back-rotation occurs is

Clast pressure
The Eshelby conjecture holds for circular as well as any kind of elliptical inclusion subjected to homogeneous far-field stresses. To verify this we derive the pressure expression within the clast as Since only constants are involved, this expression yields a constant pressure value within the inclusion. We can also see that the horizontal simple shear (ε = 0 and α = 0) does not cause a pressure perturbation within the elliptical inclusion. The same is true for a pure shear-only case where the boundary condition is 45 • inclined. In the limit of a circular inclusion,ρ c → ∞, the pressure perturbation is likewise 0. It is also obvious that in simple as well as pure shear, the rotating inclusion goes through two maxima and two minima in one full rotation. For the combined pure and simple shear an interesting parameter is the inclination angle at which, for given boundary condition amplitudes, the maximum absolute pressure deviation from the background state occurs. Taking the derivative of p c versus α, equating to 0 and solving for α we obtain For simple shear only this yields α = 45 • and for pure shear only it gives α = 0 • . Substituting eq. (68) into eq. (67) results in the actual expression for the maximum, absolute pressure perturbation within the elliptical inclusion: Deriving the maximum pressure perturbation occurring in the rigid (μ c → ∞) elliptical inclusion in pure shear, we obtain  This shows that the maximum pressure perturbation, equal to the overpressure, is roughly equal to t/2 times the characteristic far-field matrix stress value.

Pressure around rigid elliptical clast
Substituting the complex potentials of the matrix into the Joukowski transformation adjusted pressure expression eq. (62), we obtain the corresponding pressure perturbation values: This expression is again closely related to the matrix pressure expression around the circular clast (eq. 34). Namely, if we set α = 0,ε = 0 and investigate the small aspect ratio case,ρ c → ∞, the two expressions are identical. Settingρ =ρ c and α = 0 the expression for the matrix pressure around the clast-matrix interface subjected to ellipse long axis parallel far-field flow can be derived. Results for different t and µ c in simple shear are shown in Fig. 7, and pure shear in Fig. 8. In order to synchronize the plots, we scale the pressure values by the characteristic far-field matrix stress. This is µ mγ in simple shear and −2µ mε in pure shear. The minus sign takes care of the pressure sign convention.
In all four cases the maximum pressure perturbation around the circular clast is twice the matrix stress value. Inversion of the clast-matrix viscosity contrast causes a sign flip of the pressure perturbation. The behaviour of the elliptical inclusion is more complicated. In simple shear the competent elliptical clast causes progressively less pressure perturbation with increasing aspect ratio, the inverse is true for the weak inclusion. An additional effect that can be observed is that the more elongated the elliptical inclusion in simple shear, the closer the pressure extrema are to the tips.
Interestingly the effect of increasing the aspect ratio of the competent inclusion is opposite in simple shear and pure shear. In pure shear the pressure perturbation around the competent elliptical inclusion grows with increasing aspect ratio. It appears that the maximum pressure (at the tips, θ = 0) grows linearly with t. On the other hand, the maximum absolute pressure perturbation driven by the presence of a weak inclusion first increases with increasing t, but then decreases again. The viscosity contrast limits of the analytical expression for the maximum pressure perturbation generated by ellipse parallel pure shear in the tips confirm this qualitative observation:  Complete pressure field Fig. 9 illustrates examples of 2-D pressure fields in and around elliptical inclusions under ellipse-parallel pure shear conditions. Again, it can be observed (Figs 9a and b) that going from a strong to a weak inclusion leads to a flip in the pressure field. The tips of the weak elliptical inclusion become regions of relative extension, although they lie in the direction of maximum compression. The effect of increasing aspect ratios is that the elliptical inclusion focuses the pressure concentrations near the tips, although the inclination of the boundary condition does not coincide with the major axis of the ellipsoidal inclusion (Figs 9c and 7b).

Clast maximum shear stress
The complete confirmation of the Eshelby conjuncture is that with eq. (63) the following expression for the elliptical inclusion can be derived: Together with the fact that we have already shown that the pressure inside the inclusion is always homogenous, it can be deduced from eq. (73) that all stress components are homogenous since only coordinate-independent constants are involved. The corresponding expression for the C 2003 RAS, GJI, 155, 269-288 For a clast in a shear zone, where generallyμ c t, we can derive a simpler expression for the maximum τ : τ c 4µ m |γ | = (t + 1) 2 8t .
For clasts with pressure-insensitive visco-plastic rheology, eq. (76) can be used to determine whether the clast is behaving viscously or starts to yield. It can be seen that at least an aspect ratio of 6:1 is required to produce a maximum shear stress value within the clast that exceeds the far-field shear stress value (4µ m |γ |).

U N D E R S TA N D I N G C L A S T S I N S H E A R Z O N E S
Clasts in shear zones are a classic example of heterogeneities in rocks. While the classical solution for this problem has been derived long ago (Jeffery 1922), many questions remain open. For example, it is unclear how isolated, rigid clasts in a shear zone can develop a shape-preferred orientation (SPO, systematic, stable orientation of clasts), when Jeffery's theory predicts that they should keep rotating as long as the far-field flow is maintained (e.g. Pennacchioni et al. 2001). We argue (Schmid 2002;Schmid & Podladchikov 2003b) that the presence of a weak rim between the rigid clast and the matrix or slipping interface conditions is the reason that SPO develops. In order to understand the influence of a weak rim, it is necessary to know the behaviour of 'end-member' clasts: very competent or very weak, embedded in the same matrix (Figs 10a and b). As demonstrated in the previous sections an inversion of the viscosity contrast always leads to a complete reversal in the pressure perturbation. Now, if a weak rim around the circular inclusion is introduced (Fig. 10c) then the question arises as to whether the weak rim or the competent clast is dominating the effective behaviour of the combined clast-rim couple. Using the pressure perturbation in the matrix as a proxy, it is clear from Fig. 10(c) that the clast-rim couple looks to the outside world as if it were a weak inclusion, i.e. the competent clast may be ignored. We have also shown in this paper that weak elliptical inclusions have the potential to back rotate against the applied shear sense towards a stable position. If these two characteristics are combined, it becomes clear how and why isolated clasts can develop an SPO (Schmid & Podladchikov 2003a).
Another interesting point concerns the development of pressure shadows. It has always been assumed that the pressure perturbation around a clast corresponds to Fig. 10(a) and that new material develops in the quadrants of relative extension. However, at the same time it is quite likely that a weak rim of fine-grained material develops around the clast. As shown in Fig. 10(c) this leads to an inversion of the pressure field perturbation in the matrix and the development of pressure shadows should move into the other quadrants. If this is indeed the case, then pressure shadows may contain misleading shear sense information. Naturally, such a hypothesis should be backed up by other methods and data: field, experimental and numerical. Nevertheless, although only covering simple cases, the provided closed-form solutions are valuable and flexible instruments for studying the basic characteristics of heterogeneous rocks.

D I S C U S S I O N A N D C O N C L U S I O N
We have presented the complete set of possible closed-form analytical solutions for isolated, deformable elliptical inclusions subjected to general shear, far-field flows. The aim was to demonstrate how the essential parameters can be extracted from the two complex potentials that form the solution based on Muskhelishvili's method and how to apply these solutions to geological problems.
The rheologies used in this work are all Newtonian. It is clear that for geological purposes, an adequate rheology would be more complex, such as a power law. The method as presented here, is applicable to incompressible all-elastic and all-viscous problems. As a result of the established validity of Eshelby's conjecture for inclusions with perfect bonding to the matrix, it is possible to introduce more complex rheologies for the inclusion since the uniformity of the inclusion values largely facilitates the analysis. For example, it is possible to introduce a power-law or visco-elasto-plastic rheology for the inclusion. However, for the matrix, or even the layer between the inclusion and the matrix, this is not possible. Yet, it must be emphasized that the analytical solution of the Newtonian case is of utmost importance.
(1) The presented analytical solutions allow for verification, testing and tuning (e.g. necessary mesh resolution) of numerical codes for the Newtonian case. This should be a pre-requisite before investigating non-linear rheologies.
(2) Only analytical solutions allow for elegant handling of large parameter spaces. Especially, if the characteristics of the system are unexpected, as shown here, analytical solutions are the fastest possible way to capture the systematics and to train the intuition. Based on the understanding obtained from the analytical solutions, the analysis can be continued with combined asymptotic and numerical methods.
(3) Generally, power-law results are a variation of Newtonian results, preserving many of the Newtonian characteristics. Examples include mantled porphyroclasts (Schmid 2002) and folding of a Newtonian layer (Biot 1961) compared with a power-law layer (Fletcher 1974).
(4) Analytical solutions allow for direct and precise analysis of key parameters throughout the domain. With numerical models, this can only be achieved through systematic runs that may involve thousands of experiments.
Our application of the provided solutions to a geological problem examines the (mantled) clast in a shear zone. The results have important implications concerning the: (1) distribution and amplitudes of the driving forces of phase changes, metamorphic reactions and deformation mechanisms; (2) estimation of clast-matrix viscosity contrasts; (3) kinematic clast behaviour and development of shape-preferred orientations; (4) strain estimation from clasts in natural shear zones; and (5) general understanding of the clast as a function of geometry, present phases and viscosity contrasts. In particular, we have shown that: (1) the driving forces for metamorphic reactions and deformation do not depend on clast size, but on clast shape; (2) weak and strong clasts cause opposite trends of the dynamic parameters in and around the clasts; (3) nevertheless, in simple shear the rotation rate of a rigid circular clast is not affected by the matrix viscosity, nor by the presence of a strong or weak rim around the clast, hence the interpretation of kinematic and dynamic parameters from field data must be carefully separated; and (4) most analysed parameters approach the infinite viscosity contrast limits rapidly and viscosity contrasts of 10:1 between clast and matrix already cause parameter values that show less than 10 per cent deviation from the infinite viscosity contrast limits. The latter greatly affects the analysis of deformed particles in rocks for the interpretation of viscosity contrast. Small (<10) values for natural viscosity ratios, such as resulting from the study of deformed conglomerates (Treagus & Treagus 2002), may therefore be due to the limitations of the employed methods.

A C K N O W L E D G M E N T S
This research was supported by the ETH Zurich, grant TH 0-20650-99. We are grateful to Harro Schmeling for his helpful review comments. We wish to thank Jean-Pierre Burg, Boris Kaus and Stefan Schmalholz for helpful and inspiring discussions.