A Multiscale Model of Partial Melts 2: Numerical Results
Abstract
In a companion paper, equations for partially molten media were derived using twoscale homogenization theory. This approach begins with a grain scale description and then coarsens it through multiple scale expansions into a macroscopic model. One advantage of homogenization is that effective material properties, such as permeability and the shear and bulk viscosity of the twophase medium, are characterized by cell problems, boundary value problems posed on a representative microstructural cell. The solutions of these problems can be averaged to obtain macroscopic parameters that are consistent with a given microstructure. This is particularly important for estimating the “compaction length” which depends on the product of permeability and bulk viscosity and is the intrinsic length scale for viscously deformable twophase flow.
In this paper, we numerically solve ensembles of cell problems for several geometries. We begin with simple intersecting tubes as this is a one parameter family of problems with well known results for permeability. Using this data, we estimate relationships between the porosity and all of the effective parameters with curve fitting. For this problem, permeability scales as , , as expected and the bulk viscosity scales as , , which has been speculated, but never shown directly for deformable porous media. The second set of cell problems add spherical inclusions at the tube intersections. These show that the permeability is controlled by the smallest pore throats and not by the total porosity, as expected. The bulk viscosity remains inversely proportional to the porosity, and we conjecture that this quantity is insensitive to the specific microstructure. The computational machinery developed can be applied to more general geometries, such as texturally equilibrated pore shapes. However, we suspect that the qualitative behavior of our simplified models persists in these more realistic structures. In particular, our hybrid numerical–analytical model predicts that for purely mechanical coupling at the microscale, all homogenized models will have a compaction length that vanishes as porosity goes to zero. This has implications for computational models, and it suggests that these models might not resist complete compaction.
Gindraft=false \authorrunningheadSIMPSON ET AL. \titlerunningheadA MULTISCALE MODEL OF PARTIAL MELTS \authoraddrG. Simpson, Department of Mathematics, University of Toronto, Toronto, ON M5S 2E4, Canada. (simpson@math.toronto.edu) \authoraddrM. Spiegelman, Department of Applied Physics and Applied Mathematics, Columbia University, New York, NY 10027, USA. LamontDoherty Earth Observatory, Palisades, NY 10964, USA. (mspieg@ldeo.columbia.edu) \authoraddrM. I. Weinstein, Department of Applied Physics and Applied Mathematics, Columbia University, 200 Mudd, New York, NY 10027, USA. (miw2103@columbia.edu)
1 Introduction
Partially molten regions in the Earth’s asthenosphere (e.g. beneath midocean ridges or subduction zones) are usually modeled as a viscously deformable permeable media. Such models are typically composed of macroscopic equations for the conservation of mass, momentum and energy of each phase. In our companion paper, [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein], we derived several systems of governing equations for partially molten systems using homogenization. Briefly, we began with a grain scale description of two interpenetrating fluids, each satisfying the Stokes equations, coupled by their common interface. Several different macroscopic models could then be coarsened from this microscopic description, depending on our assumptions on the velocities, viscosities, and grain scale geometry.
An important feature of this approach is that macroscopic properties such as permeability, shear viscosity, and bulk viscosity naturally appear in the macroscopic equations, even if they are not defined at the grain scale. In particular, permeability and bulk viscosity are properties of the twophase aggregate, not the volume averages of small scale variations. In contrast, previous work on the magma problem, including [McKenzie(1984), Bercovici and Ricard(2003), Bercovici and Ricard(2005), Bercovici(2007), HierMajumder et al.(2006)HierMajumder, Ricard, and Bercovici, Ricard(2007)], began with models much larger than the grain scale. There, the viscosities, permeability, and other closures were assumed and justified from other results. In contrast, homogenization derives these properties selfconsistently.
While homogenization techniques appropriately inserts the constitutive relationships into the macroscopic equations, connecting them to the microstructure requires the solution of specific “cell problems.” For the physical system derived in [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein], there are actually a series of ten Stokes problems for both fluid and solid posed on the microscale (which can be reduced to four for microstructures with sufficient symmetry).
In this paper, we numerically explore the cell problems to extract parameterizations of the various effective parameters in terms of porosity, a simple measurement of the microstructure. A feature of this work is to derive constitutive relationships for the bulk viscosity as a function of porosity, which is essential for describing compactible permeable media. In particular, we show that for a range of simple pore structure the effective bulk viscosity is related to the porosity as.
Our results and some additional theory suggest that this scaling is insensitive to the specific pore geometry.
Since permeability and bulkviscosity can both be consistently related to the same microstructure, these calculations also allow us to study the “compaction length” [McKenzie(1984)]. At sufficiently low melt concentrations, it is approximately
which depends on the product of the derived permeability and bulk viscosity , both of which depend on the porosity. The compaction length is the intrinsic length scale in magma dynamics. This work suggests that under solely mechanical deformation,
implying that no mechanical mechanism prevents a region from compacting to zero porosity. This result also places strong resolution constraints on computational models of magma migration which may require a regularization for small porosities.
An outline of this work is as follows. In Section 2, we review several models of partially molten rock and highlight the constitutive relations. In Section 3, we demonstrate the process of assuming a cell geometry to extract computational closures for the macroscopic system. We revisit these closures in Section 4 for more general geometries to assess their robustness. Finally, in Section 5, we combine our numerics with the equations and examine the implications.
2 Review of Equations and Constitutive Relations
2.1 Macroscopic Equations
In [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein], we showcased three models for momentum conservation in a partially molten medium. They were distinguished by the assumed scalings for the relative velocities and viscosities between the fluid and solid phases, along with the connectivity of the pore network. One of them, dubbed BiphasicI, was given by the equations:
(1a)  
(1b)  
(1c) 
Notation for his model may be found in Table 1. In particular, and are the emergent permeability and bulk viscosity. is an auxiliary, tensorial, shear viscosity capturing grain scale anisotropy. All are related to the aforementioned cell problems. We highlight this case because (1a – 1c) is nearly identical to the models of McKenzie, Bercovici, Ricard, and others in the absence of melting and surface physics.
Symbol  Meaning 

Compaction length  
Strain rate tensor,  
Supplementary anisotropic viscosity derived by homogenization  
Porosity  
Permeability tensor of the matrix derived by homogenization  
Isotropic permeability of the matrix derived by homogenization  
Shear viscosity of the melt  
Shear viscosity of the matrix  
Macroscopic (fluid) pressure derived by homogenization  
Melt density  
Matrix density  
Mean density,  
Macroscopic fluid velocity derived by homogenization  
Macroscopic solid velocity derived by homogenization  
Bulk viscosity of the matrix derived by homogenization 
2.2 Constitutive Relations
The constitutive relations for the permeability and viscosity are fundamental to the dynamics of these models. Indeed, they are the source of much nonlinearity and it is useful to review some proposed closures. Notation for this appears in Table 2.
Symbol  Meaning 

Critical porosity for activation of viscosity  
Permeability of the matrix  
Shear viscosity of the matrix in the presence of melt  
Bulk viscosity of the matrix 
At low porosity, it is common to relate permeability, , to porosity by a power law, . Estimates of vary, , [Scheidegger(1974), Bear(1988), Dullien(1992), Turcotte and Schubert(2002), McKenzie(1984), Doyen(1988), Cheadle(1989), Martys et al.(1994)Martys, Torquato, and Bentz, Faul et al.(1994)Faul, Toomey, and Waff, Faul(1997), Faul(2000), Koponen et al.(1997)Koponen, Kataja, and Timonen, Wark and Watson(1998), Wark et al.(2003)Wark, Williams, Watson, and Price, Cheadle et al.(2004)Cheadle, Elliott, and McKenzie]. For partially molten rocks, the exponent is better constrained by both analysis and experiment to .
For the matrix shear viscosity, [Hirth and Kohlstedt(1995a), Hirth and Kohlstedt(1995b), Kohlstedt et al.(2000)Kohlstedt, Bai, Wang, and Mei, Kelemen et al.(1997)Kelemen, Hirth, Shimizu, Spiegelman, and Dick, Kohlstedt(2007)] experimentally observed a melt weakening effect, which they fit to the curve:
(2) 
is the shear viscosity of the solid matrix in the presence of melt; is the shear viscosity a meltfree matrix. In [Bercovici et al.(2001)Bercovici, Ricard, and Schubert, Bercovici and Ricard(2003)] and related works, the viscosity is weighted by , which is also present in (1a). Reiterating, is a fitting of experimental data. Regardless, the shear viscosity is taken to be isotropic, and porosity weakening in other models.
Lastly, the bulk viscosity, , is often taken as , , though is usually either zero or one. has often been used because the variation of bulk viscosity with porosity is poorly constrained by observations. [Scott and Stevenson(1984)] and others invoked the bore hole studies of ice by [Nye(1953)] to justify . In that work, Nye considered the dynamics of individual spherical and cylindrical voids in an infinite medium. [Taylor(1954), Prud’homme and Bird(1978)] computed , in the limit of small porosity, using models of incompressible fluids mixed with gas bubbles. [Schmeling(2000)] also employed , citing studies for the analogous question of the effective bulk modulus of a fluid filled poroelastic medium. Those results rely on the selfconsistent approximation methodology, treated in [Torquato(2002)]. Models in [Bercovici et al.(2001)Bercovici, Ricard, and Schubert, Bercovici and Ricard(2003)] and related works possess a term that functions as a bulk viscosity, though it has a very different origin. [McKenzie(1984)] suggested using a metallurgical model of spheres due to [Arzt et al.(1983)Arzt, Ashby, and Easterling] which considered a range of mechanisms for densification of powders under hot isostatic pressing. For viscous compaction mechanisms, the Arzt result recover a bulk viscosity proportional to for Newtonian viscosities. However, for very small porosities, they suggest that surface diffusion effects become important and imply that . Finally, [Connolly and Podladchikov(1998)] invoke an assymetric bulk viscosity that is weaker during expansion than during compaction, which is motivated by the deformation of brittle crustal materials. However, it is unclear if this rheology is relevant to hightemperature/highpressure creeping materials as are expected in the mantle. In general, it is expected that the bulkviscosity should become unbounded as the porosity reduces to zero as the system then becomes incompressible.
2.3 Cell Problems
In homogenization, a medium with fine scale features is modeled by introducing two or more spatial scales. As discussed in [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein], the direct approach makes multiple scale expansions of the dependent variables, letting them depend on both the coarse and fine length scales.
The two characteristic lengths in our model are , the macroscopic scale, and , the grain scale. Their ratio,
is the parameter in which we perform the multiple scale expansions. In the process of performing the expansions and matching powers of , we are mathematically constrained to solve a collection of auxiliary cell problems posed on the fine scale. Notation for this setup is provided in Table 3.
Symbol  Meaning 

Ratio of microscopic and macroscopic length scales,  
Strain rate tensor,  
Interface between melt and matrix within cell  
Grain length scale  
Macroscopic length scale  
Gradient taken with respect to the argument  
Divergence taken with respect to the argument  
Macroscopic region containing both melt and matrix  
Coordinate within the cell,  
The unit cell  
Portion of unit cell occupied by melt  
Portion of unit cell occupied by matrix  
Pressure of the cell problem for a unit forcing on the divergence equation 
It is analytically advantageous to approximate the mixture as a periodic medium; such a configuration appears in Figure 1. , the macroscopic region containing both matrix and melt is periodically tiled with scaled copies of the cell, scaled to unity in Figure 2. The cell problems are posed within , the matrix portion of the cell, and , the melt portion for the cell. and meet on interface .
Generically, the cell problems take the form:
(3a)  
(3b)  
(3c) 
which are generally, compressible Stokes flow problems for the microscopic solid and melt velocities and pressures . Here, is the gradient, is the divergence, and the strainrate operator defined on the fine length scale . A more complete description is given in [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein]. , , , and are prescribed forcing functions on the relevant portion of , either or . The solution, , is periodic on the portion of the boundary not intersecting the interface . The cell problems may be interpreted as the fine scale response to an applied stress on a unit cell of either the melt or the matrix. The material parameters – , , and – are then defined as the cell average of an appropriate manipulation of these pressures and velocities defined at the pore scale. In particular, an average of the solid pressure in one problem determines the bulk viscosity and certain averages of melt velocities in another problem determine the permeability. Again, we emphasize that the material parameters are not volume averages of the same parameters at the fine scale.
Symbol  Meaning 

Velocity of the cell problem for a unit shear stress forcing on the solid in the component of the stress tensor  
Unit vector in the i–th coordinate,  
Effective porosity, portion of the porosity in which there is appreciable flow. .  
Velocity of the cell problem for a unit forcing on the fluid in the direction  
First component of the velocity from the cell problem with unit forcing on the fluid in the direction  
Velocity of the cell problem for a unit forcing on the divergence equation  
Volume average of a quantity over the melt portion of a cell,  
Volume average of a quantity over the matrix portion of a cell,  
Pressure of the cell problem for a unit shear stress forcing on the solid in the component of the the stress tensor  
Pressure of the cell problem for a unit forcing on the fluid in the direction  
Pressure of the cell problem for a unit forcing on the divergence equation 
3 Effective Parameters: Intersecting Tube Geometry
To connect the effective material properties – , , and – with the porosity, we must solve the cell problems and extract a parameterization. Unfortunately, the solution of the Stokes equations in a generic three dimensional domain lacks an analytic representation. Thus we compute the solutions to the cell problems numerically and fit the results to appropriate parameterized constitutive models. Notation for this section is summarized in Table 4.
As a first example of numerically closing the constitutive relations in a homogenization based model, we study the cell domain of triply intersecting cylinders, pictured in Figure 3. The fluid occupies the cylinders while the solid is the complementary portion of the cube. While this geometry is an oversimplification of real poregeometries, it serves as proof of concept of the unified, selfconsistent homogenization algorithm. Additionally, it serves as a numerical check since the permeability for this problem is expected to scale as . In Section 4 we examine a generalization of this geometry.
In what follows, we emphasize that the information given by the parameterizations should be read with a certain skepticism. The order of magnitude and signs of the coefficients and exponents are of greater utility than the particular numbers.
We remark here that for the intersecting tube geometry, the tube radius scaled to the cellsize, , can be explicitly related to the porosity:
(4) 
We solve our problems using finite elements on unstructured meshes using the open source libraries from the FEniCS and PETSc projects [[, e.g.]]DupHof2003,KirLog2007,Log2007,petscuserref,petscwebpage, with meshes generated using CUBIT [Sandia Corporation(2008)]. Details of this method and numerical benchmarks are given in Appendix B.
3.1 Effective Permeability
The first cell problem we treat is for permeability. The equations are:
(5a)  
(5b)  
(5c) 
where is a threedimensional velocity and is a scalar pressure, for each . These are equations (64a – 64c) in [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein] and describe the motion of the melt through the porous matrix. In general, the permeability is the second order tensor:
The symmetry properties of the domain simplify this to:
Thus, it is sufficient to compute the case . Then , the permeability of the matrix, in (1b) is
(6) 
Darcy’s Law and permeability have been studied by many techniques, including homogenization; we refer the reader to the references in Section 2.2. We study it here to understand how the permeability behaves in concert with the other constitutive relations as the microstructure varies. This also serves as a benchmark problem for our software; see Appendix B.
As noted, porosity and permeability are often related by a power law, , with . To motivate such a relation, we turn to a toy model, as presented in [Turcotte and Schubert(2002)]. The melt is assumed to be in Poiseuille flow through triply intersecting cylinders. Additionally, the cylinders have small radii; it is a low porosity model. The permeability of such a system is
(7) 
Other simple models are developed in [Scheidegger(1974), Bear(1988), Dullien(1992)].
We now fit our computed permeabilities, , to porosity by such a relation. For the tube domains, the least squares fit is
(8) 
This curve and the data appear in Figure 4. The fit matches expectations of an prefactor and an exponent . The error in (8) is the associated 95% confidence interval. We report these intervals in all regressions, though they rely on the specious assumption that error in our synthetic data is normally distributed.
3.2 Effective Bulk Viscosity
The effective bulk viscosity is related to the solution of:
(9a)  
(9b)  
(9c) 
where is a threedimensional velocity and is a scalar pressure. These are equations (B1a – B1c) in [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein] and are associated with the compaction of the matrix. The effective bulk viscosity is then defined as
(10) 
The dependence of the effective bulk viscosity of partially molten rock as a function of porosity is the most poorly constrained of the material properties. This is partly due to the difficulties in constructing an experiment that will measure it as a function of porosity independently of the shear viscosity [McKenzie(1984), Kelemen et al.(1997)Kelemen, Hirth, Shimizu, Spiegelman, and Dick, Stevenson and Scott(1991)]. As mentioned in Section 2.2, a bulk viscosity has often appeared in the literature.
Two toy models for the bulk viscosity of an incompressible fluid seeded with compressible gas bubbles were formulated by [Taylor(1954), Prud’homme and Bird(1978)]. They relate the bulk viscosity to the porosity as:
(11a)  
(11b) 
Taylor’s expression, (11a), relied on a single inclusion model for a gas bubble in an infinite medium. (11b) is derived by considering a sphere of fluid with a gas filled spherical cavity, and seeking the bulk viscosity of a compressible fluid that will give rise to the same radial stress for specified boundary motion. The is due to Prud’homme and Bird restricting their model to a finite volume. This factor also appears in the proposed bulk viscosity of [Schmeling(2000)]. These expressions have motivated the use of in macroscopic models of partial melts, although they actually arise from a different problem of a compressible inclusion in an incompressible fluid, rather than the divergent flow of two interconnected incompressible fluids. Surprisingly, when we consider a simple toy problem to approximate the the full cell calculation, we find that the expressions are identical.
Appendix A provides details of this toy problem, which is related to equations (9a – 9c), and gives the solution
from which we get
This motivates trying to numerically fit to , expecting and to be close to unity. Indeed, the data, plotted in Figure 5, fits the curve
(12) 
which is quite similar to the scaling of the toy model.
3.3 Supplementary Anisotropic Viscosity
We now examine the cell problem related to the supplementary viscosity , a fourth order tensor. The equations are:
(13a)  
(13b)  
(13c) 
where is a threedimensional velocity vector and is a scalar pressure, for each pair , and . These are equations (B2a – B2c) from [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein] and are tied to tensorial surface stresses applied on the matrix. Using the solutions ,
(14) 
Because of symmetry, we need only consider two problems: and , corresponding to normal stress and shear stress in each direction.
Although we have no toy problem as motivation, proved to be satisfactory. First, we study the problem , a uniaxial stress problem. For the tubes, we fit
(15) 
This vanishes as and as and it is nearly linear at small porosity. The curves and the data are plotted in Figure 6.
There is also the simple shear stress problem, . For this, we fit
(16) 
Plots for this are given in Figure 7.
4 Generalization of the Cell Domains
Regrettably, Earth materials are not as trivial as intersecting cylinders. Even an idealized olivine grain is a tetrakaidekahedron, pictured in Figures 8. As depicted, some fraction of the melt lies along the triple junctions and some is at the quadruple junctions. Other examples of idealized, texturally equilibrated arrangements appear in [von Bargen and Waff(1986)] and [Cheadle(1989), Cheadle et al.(2004)Cheadle, Elliott, and McKenzie]. These methods are also amenable to studying random media. One could compute relations based on ensembles of randomly generated cell domains, solving all relevant cell problems on the ensemble.
Motivated by Figure 8, we explore a simple generalization of the tube geometry by adding a sphere of independent radius at the intersection, as in Figure 9. This retains the symmetry of the previous model, but adds a second parameter, allowing multiple geometries for the same porosity. The sphere captures some aspect of the pocket at the quadruple junctions. The sphere radius, , and the tube radius, , are related to the porosity by the equation:
(17) 
We shall refer to it as the sphere+tube geometry.
We now repeat the computations of Section 3 on the sphere+tube geometry to assess the sensitivity of the effective parameters to cell geometry. We also perform computations on domains where the fluid occupies an isolated sphere at the center of a cube. Though this is disconnected, it provides useful information.
4.1 Permeability Revisited
Parameterizing the porositypermeability relation for this generalization is quite challenging. In contrast to the tube geometry, the data points of the sphere+tube geometry, plotted in Figure 4, do not collapse onto a curve. There is some positive correlation between permeability and porosity, and for a given porosity, the permeability of the equivalent tube geometry is an upper bound.
To better understand the trend, we examine the computed flow fields in Figures 10 and 11. These plot the velocity magnitude on two fluid domains with the same tube size, but different sphere sizes. Most of the flow is within the tube. While there is some detrainment as it enters the sphere, the flow in the tube appears insensitive to the size of the sphere.
This motivates fitting permeability against tube radius. Indeed, an alternative to (7), is
(18) 
The tube diameter , is equivalent to , the the tube radius in the tube and sphere+tube geometries. Both data sets appear in Figure 12. This is a significant improvement over Figure 4. The least square fits are:
(19)  
(20) 
These estimates with (18); taking and scaling out , this relationship is . The data is still positively correlated with sphere radius, altering it by as as much as an order of magnitude. The deviations are greatest when both and .
That the permeability is more strongly correlated with the tube radius than the overall geometry is not surprising. [Koponen et al.(1997)Koponen, Kataja, and Timonen] discuss the notion of effective porosity, the portion of the void space where there is significant flow. Denoting our effective porosity , we seek a relation .
Given the flow fields in Figures 10 and 11 and the success with the tube radius fittings, we posit that the effective porosity is the portion of the porosity within the tubes. A twodimensional analog appears in Figure 13. Using (8), we define for the sphere+tube domains:
(21) 
[Zhu and Hirth(2003)] made a similar approximation; from [von Bargen and Waff(1986)] they construed that the permeability was controlled by the minimal crosssectional area of the pore network. A similar argument is made by [Cheadle(1989)]. In our domains, the minimal crosssectional area is .
We fit
(22) 
This appears in Figure 14. Again, deviation is highest for very large spheres with very thin tubes. Unfortunately, does not satisfy a conservation law, making it a less than ideal macroscopic quantity to track, though it does satisfy the bound .
There is still as much as an order of magnitude deviation at low porosity from relation (22). The unresolved part of the permeability for the fit is increasing in the sphere radius. This motivates trying to fit against both and another parameter. It is sufficient to fit permeability to and , resulting in
(23) 
This is plotted in Figure 15. Some deviation persists at low porosity, but it is less than an order of magnitude. Both the sphere+tube data points and the tube data points collapse onto this curve. (23) is also consistent with (8), the fit of porosity agaist permeability for the tube geometry. Taking for the tubes, (23) becomes
(24) 
This is similar to the most general permeability relationships, formulated in [Scheidegger(1974), Bear(1988)]:
(25) 
By including both and in (23), we capture some aspect of the pore shape.
4.2 Bulk and Supplementary Viscosities Revisited
In contrast to the permeability problem, the effective bulk viscosity and supplementary anisotropic viscosity are quite robust to the domain distortion. As before, we fit to , expecting and to be close to unity. For the sphere+tube and the sphere geometry, the least squares fits are:
(26a)  
(26b) 
The data and these fits are plotted appear in Figure 5. The spherical geometry appears to be an upper bound on the bulk viscosity for a given porosity. We also remark that the prefactors vary by less than an order of magnitude amongst the different domains. This is a strong endorsement of an effective bulk viscosity not only for small porosity, but also for moderate porosities .
The supplementary anisotropic viscosity terms are also robust under this geometric perturbation. For the problem , the two new geometries fit:
(27a)  
(27b) 
The curves and the data are also plotted in Figure 6. For , the spread amongst the three geometries is less than an order of magnitude.
Similar results are found in the problem. The two new additional domains are fit with:
(28a)  
(28b) 
The sphere+tube data is bounded between the sphere data and the tube data; the spread is less than an order of magnitude.
5 Discussion and Open Problems
We have successfully parameterized the macroscopic parameters on ensemble of domains for the simple pore geometries. These can now be consistently used with the macroscopic model given by equations (1a – 1c). We now review and discuss our computations, both independently of and together with the macroscopic equations.
5.1 Sensitivity to Geometry
As demonstrated by our computations in Sections 3 and 4, the effective parameters demonstrate a variety of sensitivities to the geometry. Permeability, for our cell geometries, can be bounded by porosity, , where , as seen in Figure 4. But in general, it cannot be expressed as a function of a single parameter, such as porosity. However, if the shapes are constrained by textural equilibration, as in [Cheadle(1989)], there will be a single variable parameterization for each dihedral angle. Solving the cell problems on these shapes, and assessing their sensitivity to the dihedral angle, is an important open problem.
In contrast, the bulk viscosity appears to be rather insensitive to the geometry, scaling as . We see this from the variation, or lack thereof, between the data and fits for the isolated spheres and the triply intersecting cylinders, in Figure 5. Though these are entirely different geometric structures, there is less than an order of magnitude of variation in the computed . These too merit computation on the texturally equilibrated shapes. Randomly generated geometries may also be of interest.
Both components of the supplementary anisotropic shear viscosity appear to be insensitive to geometry, with less than an order of magnitude amongst the three geometries. However, because the problems are driven by surface stresses, it may be that a more anisotropic shape could alter these scalings. Studying them on randomly generated shapes may provide insight on the role of grain scale anisotropy.
5.2 Bulk Viscosity
Perhaps our most significant result is the selfconsistent bulk viscosity, arising from a purely mechanical model of partially molten rock. A spatially varying bulk viscosity is quite important. Indeed, significant differences in dynamics were noted between the solutions of the [McKenzie(1984)] model and the model in [Ricard et al.(2001)Ricard, Bercovici, and Schubert]. The authors point to the use of a constant bulk viscosity in the McKenzie model as the source of the discrepancy.
Prior to [Ricard et al.(2001)Ricard, Bercovici, and Schubert], [Schmeling(2000)] remarked that the dependence has an important impact on the compaction length. For , the bulk viscosity is two orders of magnitude greater than the shear viscosity. While many studies took and to be the same order, this higher bulk viscosity leads to a compaction length an order of magnitude greater. [Ricard(2007)] made a similar observation on the impact of variable bulk viscosity on the compaction length.
Schmeling also commented that this variable bulk viscosity could induce melt focusing towards the axis in his plume simulations. This is an additional nonlinearity that may be important to geophysical problems. Many studies relying on the McKenzie model employed a constant bulk viscosity, including [Richter and McKenzie(1984), Spiegelman and McKenzie(1987), Spiegelman(1993a), Spiegelman(1993b), Spiegelman(1993c), Aharonov et al.(1997)Aharonov, Spiegelman, and Kelemen, Spiegelman(1996), Kelemen et al.(1997)Kelemen, Hirth, Shimizu, Spiegelman, and Dick, Spiegelman et al.(2001)Spiegelman, Kelemen, and Aharonov, Katz et al.(2004)Katz, Spiegelman, and Carbotte, Spiegelman et al.(2007)Spiegelman, Katz, and Simpson]. [Spiegelman and Kelemen(2003), Spiegelman(2003)] used a bulk viscosity, with with , the exponent in the permeability relationship to prevent the system from compacting to zero between the reactive channels. It would be interesting to revisit these problems with a bulk viscosity.
5.3 Compaction Length
We now combine (1a – 1c), our leading order equations derived by multiple scale expansions, with our computed constitutive relations. For our numerical estimates, (23), (26a–26b), (15–27b), (27a), and (16–28a) are approximately:
is a constant, is an constant, and is a constant fourth order tensor. Under these assumptions, the equations for the BiphasicI model, (1a – 1c), simplify:
(29)  
(30)  
(31) 
If we were to use these equations and numerically derived constitutive relations to solve a boundary value problem, the compaction length would again appear as an important length scale,
If , then and ,
(32) 
Since , we have an upper bound on the compaction length,
(33) 
Hence,
We believe (33), which constrains the compaction length by the porosity raised to a small positive power, is relatively insensitive to the geometric configuration. This follows from our alleged robustness of our effective bulk viscosity, , and the broad agreement in the porosity–permeability relationship, with .
A compaction length that vanishes with porosity has interesting consequences. For example, this compaction length scaling does not rule out the possibility that a partially molten rock could expel all fluid by mechanical means. It also does not permit the infiltration of fluid into a dry region. Understanding how any of these systems of equations transition between a partially molten region and a dry region is an outstanding question. We also note that though our scaling relationship does not forbid compaction, it does not imply it either. There may be a dynamic response that prevents the matrix from mechanically compacting to zero. Such effects were mathematically proven to exist in a onedimensional simplification of the model, without melting or freezing, in [Simpson et al.(2007)Simpson, Spiegelman, and Weinstein, Simpson and Weinstein(2008), Simpson et al.(2008c)Simpson, Weinstein, and Rosenau].
If, instead, we had concluded , with , then the compaction length would become unbounded as the melt vanished. Hence the region of deformation in the matrix needed to segregate additional fluid would also become infinite, precluding further segregation solely by mechanical processes.
Because we have parameterized the permeability with , we can explicitly see the response of the compaction length as the melt network becomes disconnected. measures the volume fraction where melt flows. As the channels close up and the melt becomes trapped and . Taking this limit in (32),
The compaction length can vanish, even if the melt fraction remains bounded away from zero. A similar conclusion could be drawn for the compaction length of McKenzie,
(34) 
Letting , if we interpret the loss of connectivity as , vanishes with nonzero porosity.
5.4 Other Physics
There are several results that were not realized by our model, and these merit discussion. We did not recover the result of [Arzt et al.(1983)Arzt, Ashby, and Easterling]. This arises in the limit of sufficiently low porosity that the primary transport mechanism is by grain boundary diffusion. Since we did not include any surface physics in our fine scale model, we should not have expected to upscale their effects. If good grain scale descriptions of these processes could be formulated, it might be possible to coarsen them via homogenization, perhaps recovering this macroscopic relation. However, as becomes unbounded more slowly than , this result does not change the basic argument about the compaction length vanishing with zero porosity.
Another relationship not captured by either our work is the experimental fit for matrix shear viscosity in the presence of melt from [Hirth and Kohlstedt(1995a), Hirth and Kohlstedt(1995b), Kelemen et al.(1997)Kelemen, Hirth, Shimizu, Spiegelman, and Dick, Kohlstedt et al.(2000)Kohlstedt, Bai, Wang, and Mei, Kohlstedt(2007)],
Our model possesses a porosity weakening mechanism; all of the viscosity terms are . However, the anisotropic part, , is not sign definite and is small compared to the isotropic component. Furthermore, there does not appear to be an exponential relation. [Hirth and Kohlstedt(1995a)] hypothesized that the presence of melt enhances grain boundary diffusion, providing a fast path for deformation through the melt. As with the bulk viscosity, this is a surface physics phenomenon not captured by our Stokes models.
5.5 Open Problems
As discussed, two important open problems are the computation of the cell problems on more realistic and general cell domains and the inclusion of surface physics into the model. The former problem is rather straightforward, requiring good computational tools for generating the domains and solving the Stokes equations on them. The latter problem is more challenging, requiring fine scale equations for these processes.
The models could also be augmented by giving the matrix a nonlinear rheology, leading to nonlinear cell problems. Though these could be solved and studied numerically, this is much more difficult as the different forcing components no longer decouple. Rather than being able to split the matrix cell problems into a bulk viscosity problem, and two surface stress problems, they would have to be done simultaneously. However, the derived parameterizations for the effective viscosities would be important to magma migration. A nonlinear matrix rheology is expected at large strain rates and was needed to computationally model physical experiments for shear bands in [Katz et al.(2006)Katz, Spiegelman, and Holtzman].
Another important physical rheology is viscoelasticity. [Connolly and Podladchikov(1998), Vasilyev et al.(1998)Vasilyev, Podladchikov, and Yuen] extended the earlier, purely viscous, models to include elastic effects. This is an important regime since it permits both short and long time scales, as is found at the asthenospherelithosphere boundary.
Appendix A Spherical Model
Here we develop a toy model for equations (9a – 9c),
whose solution yields the effective bulk viscosity,
Consider a fluid domain, , occupying a small isolated sphere at the center of the unit cube; is the complementary region. Smoothing out the exterior boundary of deforms it into a sphere. We solve the dilation stress problem on this domain. To avoid confusion, let
(35)  
(36) 
Our toy problem is posed on .
Although periodicity is no longer a meaningful boundary condition, it can be shown that the normal velocity on the periodic part of vanishes. We set the normal velocity to zero on the exterior boundary of . On the interior shell, the stress free condition remains. The equations are:
(37a)  
(37b)  
(37c)  
(37d) 
is the radius of the interior sphere. Decomposing the velocity into incompressible, , and compressible, , components, the compressible part solves:
(38) 
Let the boundary conditions on the potential be:
(39) 
Since the problem is spherically symmetric, its solution is
(40) 
The incompressible velocity must be divergence free. Again, by spherical symmetry,
To satisfy the boundary condition at , . Therefore, and . The pressure then solves , so it is constant,
Applying boundary condition (37c),
In this geometry, , and the cell averaged pressure is
(41) 
Therefore,
(42) 
(42) uses definition (10), but the ’s in each solve different problems. (42) is specific to the geometry of a sphere with cavity, while (10) is for a generic geometry occupying some fraction of the unit cube.
(41) is plotted along with data for the numerical solutions of the cell problem posed on in Figure 16. There is good agreement between (41) and these computations for porosity suggesting our deformation was reasonable. We reiterate that though our result is quite similar to that of (11a) and (11b), the origin of the underlying Stokes problem is quite different.
Appendix B Computational Methods and Results
b.1 Cell Problem Boundary Conditions
Though we could solve the cell problems as stated, with the specified boundary conditions on interface and periodic on the rest of the domain, we use the symmetry of our model problems to reduce the computational cost by a factor of eight. The symmetry properties of the solutions allow us to formulate the appropriate boundary conditions on the portion of the boundary that is not . The symmetry properties of the cell problems are summarized in Table 5. We specify these Dirichlet boundary conditions on the velocity together with the original boundary condition on . When forming the weak form of the problem, we use neutral boundary conditions, on the part of the boundary that is not . For example, in solving the permeability cell problem of Section 3.1 for and , the weak form is
(43) 
where and are test functions. The Dirichlet boundary conditions are indicated in Figure 17. On the part of the boundary that is the interface, , we have applied the no  slip condition.
Cell Problem Velocity  

b.2 Solver Algorithms
We discretize the Stokes equations for the cell problems using the P2P1 formulation described in [Elman et al.(2005)Elman, Silvester, and Wathen]. The FEniCS libraries are used to generate code for the weak forms of the equations and assemble the associated matrices and vectors, [Dupont et al.(2003)Dupont, Hoffman, Johnson, Kirby, Larson, Logg, and Scott, Kirby and Logg(2006), Kirby and Logg(2007), Logg(2007)]. These vectors and matrices are passed to PETSc and solved using algebraic Multigrid preconditioned GMRES, [Balay et al.(1997)Balay, Gropp, McInnes, and Smith, Balay et al.(2004)Balay, Buschelman, Eijkhout, Gropp, Kaushik, Knepley, McInnes, Smith, and Zhang, Balay et al.(2001)Balay, Buschelman, Gropp, Kaushik, Knepley, McInnes, Smith, and Zhang]. Domains and meshes were created with CUBIT [Sandia Corporation(2008)]. The versions of the software we used are summarized in Table 6.
Package  Version 

CUBIT  11.0 
DOLFIN(FEniCS)  0.7.2 
FFC(FEniCS)  0.4.4 
FIAT(FEniCS)  0.3.4 
HYPRE  2.0.0 
PETSc  2.3.3 
UFC(FEniCS)  1.1 
UMFPACK  4.3 
To study problems with elements and unknowns, we rely on a Stokes preconditioner employing the pressure mass matrix of [Elman et al.(2005)Elman, Silvester, and Wathen]. The Stokes system is
where is matrix corresponding to the weak form of the vector Laplacian, is the matrix corresponding to the weak form of the gradient, and is the matrix corresponding to the weak form of the divergence. This is preconditioned with an approximate inverse of
is the pressuremass matrix.
As our meshes are unstructured, the HYPRE library is used for algebraic multigrid preconditioning. In particular, we use BoomerAMG. We apply this on all of , although we could have only used this on the block, and relied on Jacobi or another light weight preconditioner for the block.
b.3 Examples and Benchmarks
As a test, we solve the permeability cell problem of Section 3.1, with the symmetry reductions, for flow past a sphere of radius . It is meshed with a characteristic size of , consisting of 119317 tetrahedrons. The results are summarized in Table 7
Method  

BoomerAMG on + GMRES  0.0447051 
JT05  0.045803 
COMSOL  0.044497 
For comparison, [Jung and Torquato(2005)] ran a time dependent problem to steady state and used an immersed boundary method with finite volumes. In our COMSOL computation, we used “fine” meshing, with 29649 elements, 134260 degrees of freedom, and a relative tolerance of 1e10 in the solver.
The objective function, , converges as we refine our mesh; see Table 8 for a comparison of different meshes for this problem.
Mesh Size  No. Cells  No. d.o.f.  

0.25  56  430  0.0492875  – 
0.125  534  2950  0.0472207  0.0419336 
0.0625  3673  17988  0.0452142  0.042492 
0.03125  26147  119317  0.0447051  0.0112597 
auto  14803  70525  0.0445419  – 
Table 9 summarizes the convergence results for flow around a sphere of radius . Again, our method appears to be quite effective.
Mesh Size  No. Cells  No. d.o.f.  

0.25  61  475  0.00763404  – 
0.125  384  2302  0.00651896  0.146067 
0.0625  2620  13370  0.00626809  0.0384831 
0.03125  19916  93011  0.00617889  0.0142308 
auto  23776  112620  0.00616139  – 
JT05  –  –  0.0064803  – 
COMSOL  –  –  0.006153  – 
As another example, we solve the dilation stress cell problem from Section 3.2 . Solved on the domain complementary to a spherical inclusion of radius , the convergence results are summarized in Table 10.
Mesh Size  No. Cells  No. d.o.f.  

0.25  73  541  46.6432  – 
0.125  514  2862  44.1747  0.052923 
0.0625  3813  18575  40.7177  0.0782575 
0.03125  29063  132115  39.4805  0.0303848 
auto  13725  64205  39.1558  – 
COMSOL  –  –  39.117074  – 
The data in Tables 8– 10 were computed with default PETSc KSP tolerances. The automatically generated mesh was constructed with the CUBIT command:
volume 6 sizing function type skeleton scale 3 time_accuracy_level 2 min_size auto max_size 0.2 max_gradient 1.3
While these convergence results are encouraging, our data is imperfect. Continuing with the dilation stress example, consider the data in Figure 18. Comparing Figures (a), (c), and (e), it would appear that the domains with smaller fluid inclusions have less well resolved pressure fields. They could likely be resolved with additional resolution. However, we use this data and believe it to be valid for several reasons:

It is preferable to have all domains meshed with the same algorithm.

While the pressure fields may not be resolved, the error appears at the interface, and we are interested in the cell average. Moreover, the relative variations about the cell average are small.

The corresponding velocity fields, with magnitudes pictured in Figures (b), (d), and (f), appear to be smooth, suggesting we are converging towards the analytical solution.

The cell averages are consistent with the trends from the better resolved cases.
b.4 Cell Problem Data
All meshes were generated using CUBIT with the command:
volume 6 sizing function type skeleton scale 3 time_accuracy_level 2 min_size auto max_size 0.2 max_gradient 1.3
Problems were solved in PETSc with a relative tolerance of and and absolute tolerance of .
Acknowledgements.
Both this paper and [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein] are based on the thesis of G. Simpson, [Simpson(2008)], completed in partial fulfillment of the requirements for the degree of doctor of philosophy at Columbia University. The authors wish to thank D. Bercovici and R. Kohn for their helpful comments. This work was funded in part by the US National Science Foundation (NSF) Collaboration in Mathematical Geosciences (CMG), Division of Mathematical Sciences (DMS), Grant DMS–05–30853, the NSF Integrative Graduate Education and Research Traineeship (IGERT) Grant DGE–02–21041, NSF Grants DMS–04–12305 and DMS–07–07850.References
 [Aharonov et al.(1997)Aharonov, Spiegelman, and Kelemen] Aharonov, E., M. Spiegelman, and P. Kelemen, Threedimensional flow and reaction in porous media: Implications for the Earth’s mantle and sedimentary basins, Journal of Geophysical Research, 102(7), 14,821–14,834, 1997.
 [Arzt et al.(1983)Arzt, Ashby, and Easterling] Arzt, E., M. Ashby, and K. Easterling, Practical Applications of Hot Isostatic Pressing Diagrams: Four Case Studies, Metall. Trans A, 14(2), 211–221, 1983.
 [Balay et al.(1997)Balay, Gropp, McInnes, and Smith] Balay, S., W. D. Gropp, L. C. McInnes, and B. F. Smith, Efficient management of parallelism in object oriented numerical software libraries, in Modern Software Tools in Scientific Computing, edited by E. Arge, A. M. Bruaset, and H. P. Langtangen, pp. 163–202, Birkhäuser Press, 1997.
 [Balay et al.(2001)Balay, Buschelman, Gropp, Kaushik, Knepley, McInnes, Smith, and Zhang] Balay, S., K. Buschelman, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, PETSc Web page, http://www.mcs.anl.gov/petsc, 2001.
 [Balay et al.(2004)Balay, Buschelman, Eijkhout, Gropp, Kaushik, Knepley, McInnes, Smith, and Zhang] Balay, S., K. Buschelman, V. Eijkhout, W. D. Gropp, D. Kaushik, M. G. Knepley, L. C. McInnes, B. F. Smith, and H. Zhang, PETSc users manual, Tech. Rep. ANL95/11  Revision 2.1.5, Argonne National Laboratory, 2004.
 [Bear(1988)] Bear, J., Dynamics of Fluids in Porous Media, Courier Dover Publications, 1988.
 [Bercovici(2007)] Bercovici, D., Mantle dyanmics past, present, and future: An introduction and overview, in Treatise on Geophysics, vol. 7, edited by G. Schubert, Elsevier, 2007.
 [Bercovici and Ricard(2003)] Bercovici, D., and Y. Ricard, Energetics of a twophase model of lithospheric damage, shear localization and plateboundary formation, Geophysical Journal International, 152(3), 581–596, 2003.
 [Bercovici and Ricard(2005)] Bercovici, D., and Y. Ricard, Tectonic plate generation and twophase damage: Void growth versus grain size reduction, Journal of Geophysical Research, 110(B3), 2005.
 [Bercovici et al.(2001)Bercovici, Ricard, and Schubert] Bercovici, D., Y. Ricard, and G. Schubert, A twophase model for compaction and damage, 1: General theory, Journal of Geophysical Research, 106(B5), 8887–8906, 2001.
 [Cheadle(1989)] Cheadle, M., Properties of texturally equilibrated twophase aggregates, Ph.D. thesis, University of Cambridge, 1989.
 [Cheadle et al.(2004)Cheadle, Elliott, and McKenzie] Cheadle, M., M. Elliott, and D. McKenzie, Percolation threshold and permeability of crystallizing igneous rocks: The importance of textural equilibrium, Geology, 32(9), 757–760, 2004.
 [Connolly and Podladchikov(1998)] Connolly, J., and Y. Podladchikov, Compaction driven fluid flow in viscoelastic rock, Geodinamica ActaRevue de Geologie Dynamique et de Geographie Physique, 11(2), 55–84, 1998.
 [Doyen(1988)] Doyen, P., Permeability, conductivity, and pore geometry of sandstone, Journal of Geophysical Research, 93(B7), 7729–7740, 1988.
 [Dullien(1992)] Dullien, F., Porous Media: Fluid Transport and Pore Structure, second edition ed., Academic Press, 1992.
 [Dupont et al.(2003)Dupont, Hoffman, Johnson, Kirby, Larson, Logg, and Scott] Dupont, T., J. Hoffman, C. Johnson, R. C. Kirby, M. G. Larson, A. Logg, and L. R. Scott, The FEniCS project, Tech. Rep. 2003–21, Chalmers Finite Element Center Preprint Series, 2003.
 [Elman et al.(2005)Elman, Silvester, and Wathen] Elman, H., D. Silvester, and A. Wathen, Finite Elements and Fast Iterative Solvers: With Applications in Incompressible Fluid Dynamics, Oxford University Press, USA, 2005.
 [Faul(1997)] Faul, U., Permeability of partially molten upper mantle rocks from experiments and percolation theory, J. Geophys. Res, 102, 10,299–10,311, 1997.
 [Faul(2000)] Faul, U., Constraints on the melt distribution in anisotropic polycrystalline aggegates undergoing grain growth, in Physics and Chemistry of Partially Molten Rocks, edited by N. Bagdassarov, D. Laporte, and A. Thompson, pp. 141–178, Kluwer Academic, 2000.
 [Faul et al.(1994)Faul, Toomey, and Waff] Faul, U., D. Toomey, and H. Waff, Intergranular basaltic melt is distributed in thin, elongated inclusions, Geophysical Research Letters, 21(1), 29–32, 1994.
 [HierMajumder et al.(2006)HierMajumder, Ricard, and Bercovici] HierMajumder, S., Y. Ricard, and D. Bercovici, Role of grain boundaries in magma migration and storage, Earth and Planetary Science Letters, 248(34), 735–749, 2006.
 [Hirth and Kohlstedt(1995a)] Hirth, G., and D. Kohlstedt, Experimental constraints on the dynamics of the partially molten upper mantle: Deformation in the diffusion creep regime, Journal of Geophysical Research, 100(B2), 1981–2001, 1995a.
 [Hirth and Kohlstedt(1995b)] Hirth, G., and D. Kohlstedt, Experimental constraints on the dynamics of the partially molten upper mantle 2: Deformation in the dislocation creep regime, Journal of Geophysical Research, 100(B2), 15,441–15,449, 1995b.
 [Jung and Torquato(2005)] Jung, Y., and S. Torquato, Fluid permeabilities of triply periodic minimal surfaces, Physical Review E, 72(5), 56,319, 2005.
 [Katz et al.(2004)Katz, Spiegelman, and Carbotte] Katz, R., M. Spiegelman, and S. Carbotte, Ridge migration, asthenospheric flow and the origin of magmatic segmentation in the global midocean ridge system, Geophys. Res. Letts, 31, 2004.
 [Katz et al.(2006)Katz, Spiegelman, and Holtzman] Katz, R., M. Spiegelman, and B. Holtzman, The dynamics of melt and shear localization in partially molten aggregates., Nature, 442(7103), 676–9, 2006.
 [Kelemen et al.(1997)Kelemen, Hirth, Shimizu, Spiegelman, and Dick] Kelemen, P., G. Hirth, N. Shimizu, M. Spiegelman, and H. Dick, A review of melt migration processes in the adiabatically upwelling mantle beneath oceanic spreading ridges, Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences, 355(1723), 283–318, 1997.
 [Kirby and Logg(2006)] Kirby, R. C., and A. Logg, A compiler for variational forms, ACM Transactions on Mathematical Software, 32(3), 417–444, 2006.
 [Kirby and Logg(2007)] Kirby, R. C., and A. Logg, Efficient compilation of a class of variational forms, ACM Transactions on Mathematical Software, 33(3), 2007.
 [Kohlstedt(2007)] Kohlstedt, D., Properties of rocks and minerals – constitutive equations, rheological behvaior, and viscosity of rocks, in Treatise on Geophysics, vol. 2, edited by G. Schubert, Elsevier, 2007.
 [Kohlstedt et al.(2000)Kohlstedt, Bai, Wang, and Mei] Kohlstedt, D., Q. Bai, Z. Wang, and S. Mei, Rheology of partially molten rocks, in Physics and Chemistry of Partially Molten Rocks, edited by N. Bagdassarov, D. Laporte, and A. Thompson, pp. 3–28, Kluwer Academic, 2000.
 [Koponen et al.(1997)Koponen, Kataja, and Timonen] Koponen, A., M. Kataja, and J. Timonen, Permeability and effective porosity of porous media, Physical Review E, 56(3), 3319–3325, 1997.
 [Logg(2007)] Logg, A., Automating the finite element method, Arch. Comput. Methods Eng., 14(2), 93–138, 2007.
 [Martys et al.(1994)Martys, Torquato, and Bentz] Martys, N., S. Torquato, and D. Bentz, Universal scaling of fluid permeability for sphere packings, Physical Review E, 50(1), 403–408, 1994.
 [McKenzie(1984)] McKenzie, D., The generation and compaction of partially molten rock, Journal of Petrology, 25(3), 713–765, 1984.
 [Nye(1953)] Nye, J., The flow law of ice from measurements in glacier tunnels, laboratory experiments and the Jungfraufirn borehole experiment, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, pp. 477–489, 1953.
 [Prud’homme and Bird(1978)] Prud’homme, R., and R. Bird, The dilational properties of suspensions of gas bubbles in incompressible newtonian and nonnewtonian fluids, Journal of NonNewtonian Fluid Mechanics, 3(3), 261–279, 1978.
 [Ricard(2007)] Ricard, Y., Physics of mantle convection, in Treatise on Geophysics, vol. 7, edited by G. Schubert, Elsevier, 2007.
 [Ricard et al.(2001)Ricard, Bercovici, and Schubert] Ricard, Y., D. Bercovici, and G. Schubert, A twophase model of compaction and damage, 2: Applications to compaction, deformation, and the role of interfacial surface tension, Journal of Geophysical Research, 106, 8907–8924, 2001.
 [Richter and McKenzie(1984)] Richter, F., and D. McKenzie, Dynamical models for melt segregation from a deformable matrix, Journal of Geology, 92(6), 729–740, 1984.
 [Sandia Corporation(2008)] Sandia Corporation, CUBIT, http://cubit.sandia.gov/, 2008.
 [Scheidegger(1974)] Scheidegger, A., The physics of flow through porous media: Third edition, University of Toronto Press, 1974.
 [Schmeling(2000)] Schmeling, H., Partial melting and melt segregation in a convecting mantle, in Physics and Chemistry of Partially Molten Rocks, edited by N. Bagdassarov, D. Laporte, and A. Thompson, pp. 141–178, Kluwer Academic, 2000.
 [Scott and Stevenson(1984)] Scott, D., and D. Stevenson, Magma solitons, Geophysical Research Letters, 11(11), 1161–1161, 1984.
 [Simpson(2008)] Simpson, G., The mathematics of magma migration, Ph.D. thesis, Columbia University, 2008.
 [Simpson and Weinstein(2008)] Simpson, G., and M. Weinstein, Asymptotic stability of ascending solitary magma waves, SIAM J. Math. Anal., 40, 1337–1391, 2008.
 [Simpson et al.(2007)Simpson, Spiegelman, and Weinstein] Simpson, G., M. Spiegelman, and M. Weinstein, Degenerate dispersive equations arising in the stuyd of magma dynamics, Nonlinearity, 20, 21–49, 2007.
 [Simpson et al.(2008a)Simpson, Spiegelman, and Weinstein] Simpson, G., M. Spiegelman, and M. Weinstein, A multiscale model of partial melts 1: Effective equations, submitted to Journal of Geophysical Research, 2008a.
 [Simpson et al.(2008b)Simpson, Spiegelman, and Weinstein] Simpson, G., M. Spiegelman, and M. Weinstein, A multiscale model of partial melts 2: Numerical results, submitted to Journal of Geophysical Research, 2008b.
 [Simpson et al.(2008c)Simpson, Weinstein, and Rosenau] Simpson, G., M. Weinstein, and P. Rosenau, On a hamiltonian pde arising in magma dynamics, DCDSB, 10, 903–924, 2008c.
 [Spiegelman(1993a)] Spiegelman, M., Flow in deformable porous media. part 1: Simple analysis, Journal of Fluid Mechanics, 247, 17–38, 1993a.
 [Spiegelman(1993b)] Spiegelman, M., Flow in deformable porous media. part 2: Numerical analysis, Journal of Fluid Mechanics, 247, 39–63, 1993b.
 [Spiegelman(1993c)] Spiegelman, M., Physics of Melt Extraction: Theory, Implications and Applications, Philosophical Transactions: Physical Sciences and Engineering,, 342(1663), 23–41, 1993c.
 [Spiegelman(1996)] Spiegelman, M., Geochemical consequences of melt transport in 2D: The sensitivity of trace elements to mantle dynamics, Earth and Planetary Science Letters, 139(12), 115–132, 1996.
 [Spiegelman(2003)] Spiegelman, M., Linear analysis of melt band formation by simple shear, Geochemistry, Geophysics, Geosystems, 4(9), 8615, 2003.
 [Spiegelman and Kelemen(2003)] Spiegelman, M., and P. Kelemen, Extreme chemical variability as a consequence of channelized melt transport, Geochemistry, Geophysics, Geosystems, 4(7), 1055, 2003.
 [Spiegelman and McKenzie(1987)] Spiegelman, M., and D. McKenzie, Simple 2D models for melt extraction at midocean ridges and island arcs, Earth and Planetary Science Letters, 83(14), 137–152, 1987.
 [Spiegelman et al.(2001)Spiegelman, Kelemen, and Aharonov] Spiegelman, M., P. Kelemen, and E. Aharonov, Causes and consequences of flow organization during melt transport: The reaction infiltration instability in compactible media, Journal of Geophysical Research, 106(B2), 2061–2077, 2001.
 [Spiegelman et al.(2007)Spiegelman, Katz, and Simpson] Spiegelman, M., R. Katz, and G. Simpson, An Introduction and Tutorial to the “McKenize Equations” for magma migration, http://www.geodynamics.org/cig/workinggroups/magma/workarea/benchmark/M%cKenzieIntroBenchmarks.pdf, 2007.
 [Stevenson and Scott(1991)] Stevenson, D., and D. Scott, Mechanics of FluidRock Systems, Annual Review of Fluid Mechanics, 23(1), 305–339, 1991.
 [Taylor(1954)] Taylor, G., The two coefficients of viscosity for an incompressible fluid containing air bubbles, Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 226(1164), 34–37, 1954.
 [Torquato(2002)] Torquato, S., Random Heterogeneous Materials: Microstructure and Macroscopic Properties, Springer, 2002.
 [Turcotte and Schubert(2002)] Turcotte, D., and G. Schubert, Geodynamics, Second ed., Cambridge University Press, 2002.
 [Vasilyev et al.(1998)Vasilyev, Podladchikov, and Yuen] Vasilyev, O., Y. Podladchikov, and D. Yuen, Modeling of compaction driven flow in poroviscoelastic medium using adaptive wavelet collocation method, Geophysical Research Letters, 25(17), 3239–3242, 1998.
 [von Bargen and Waff(1986)] von Bargen, N., and H. Waff, Permeabilities, interfacial areas and curvatures of partially molten systems: Results of numerical computations of equilibrium microstructures, Journal of Geophysical Research, 91(B9), 9261–9276, 1986.
 [Wark and Watson(1998)] Wark, D., and E. Watson, Grainscale permeabilities of texturally equilibrated, monomineralic rocks, Earth and Planetary Science Letters, 164(34), 591–605, 1998.
 [Wark et al.(2003)Wark, Williams, Watson, and Price] Wark, D., C. Williams, E. Watson, and J. Price, Reassessment of pore shapes in microstructurally equilibrated rocks, with implications for permeability of the upper mantle, Journal of Geophysical Research, 108(2050), 2003.
 [Zhu and Hirth(2003)] Zhu, W., and G. Hirth, A network model for permeability in partially molten rocks, Earth and Planetary Science Letters, 212(34), 407–416, 2003.