Textural equilibrium melt geometries around tetrakaidecahedral grains

In textural equilibrium, partially molten materials minimize the total surface energy bound up in grain boundaries and grain–melt interfaces. Here, numerical calculations of such textural equilibrium geometries are presented for a space-filling tessellation of grains with a tetrakaidecahedral (truncated octahedral) unit cell. Two parameters determine the nature of the geometries: the porosity and the dihedral angle. A variety of distinct melt topologies occur for different combinations of these two parameters, and the boundaries between different topologies have been determined. For small dihedral angles, wetting of grain boundaries occurs once the porosity has exceeded 11%. An exhaustive account is given of the main properties of the geometries: their energy, pressure, mean curvature, contiguity and areas on cross sections and faces. Their effective permeabilities have been calculated, and demonstrate a transition between a quadratic variation with porosity at low porosities to a cubic variation at high porosities.

In textural equilibrium, partially molten materials minimize the total surface energy bound up in grain boundaries and grain-melt interfaces. Here, numerical calculations of such textural equilibrium geometries are presented for a space-filling tessellation of grains with a tetrakaidecahedral (truncated octahedral) unit cell. Two parameters determine the nature of the geometries: the porosity and the dihedral angle. A variety of distinct melt topologies occur for different combinations of these two parameters, and the boundaries between different topologies have been determined. For small dihedral angles, wetting of grain boundaries occurs once the porosity has exceeded 11%. An exhaustive account is given of the main properties of the geometries: their energy, pressure, mean curvature, contiguity and areas on cross sections and faces. Their effective permeabilities have been calculated, and demonstrate a transition between a quadratic variation with porosity at low porosities to a cubic variation at high porosities.

Introduction
The physical properties of partially molten materials depend crucially on the geometry of melt at the scale of individual grains. Properties like permeability or electrical conductivity can be radically different depending on whether melt forms a connected network or not. The aim of this contribution is to better understand the controls on the geometry of melt networks in order to ultimately better understand the physical properties of partially molten materials.

The model
The model geometry consists of an infinite tessellation of tetrakaidecahedral unit cells as depicted in figure 1. In the absence of melt, a single grain occupies the whole of a unit cell, and the grain takes the shape of a plane-faced tetrakaidecahedron (truncated octahedron). The faces of the tetrakaidecahedron form the grain boundaries (solid-solid contacts). When melt is present, it is assumed that grain boundaries continue to lie along the faces of the tetrakaidecahedral cell, although the area and shape of these contacts are allowed to vary. Textural equilibrium involves the minimization of surface energy E [1,4], where γ ss and γ sl are the surface energies per unit of area of grain-grain (solid-solid) contacts and grain-melt (solid-liquid) contacts, respectively. A ss and A sl are the corresponding areas of solid-solid and solid-liquid contacts per unit cell. In this work, it will be assumed that γ ss and γ sl are isotropic and constant. The factor of 1   are on the faces of the unit cell, so per unit cell each solid-solid contact only counts for half in the total surface energy [4].
Minimization of E is performed subject to constraints, which reflect assumptions about the geometry. Here, grain centres are assumed to reside on a body-centred cubic (bcc) lattice, all grain boundaries are assumed to be planar, and all grains are identical. The unit cell shown in figure 1 is the Wigner-Seitz cell of the bcc lattice of grain centres. As a consequence of the cubic symmetry, the basic computational domain can be reduced to a single region that is 1/48th of a grain or 1/8th of a quadruple junction (figure 2).
Grain centres are fixed in these calculations, and thus there is a constant distance d between opposing square faces in the unit tetrakaidecahedral cell. The volume of the unit cell is V cell = d 3 /2 and it has area A cell = 3 4 (1 + 2 √ 3)d 2 . Each edge has length a = d/ (2 √ 2). The volume fraction of melt is prescribed, by introducing a Lagrange multiplier λ and minimizing the functional J, where V l is the volume of liquid in the domain and V * l is the target (prescribed) value. The variational problem is discretized by representing the surface by a finite-element mesh of triangles, using the SURFACE EVOLVER software [15][16][17]. Numerical optimization is used to find the unknown mesh node coordinates that minimize the surface energy subject to the given constraints. Further details on the numerical methods can be found in appendix A.
Here, H is the mean curvature of the solid-liquid surface, defined by   Thus, in textural equilibrium, the mean curvature of the solid-liquid surface is a constant. Moreover, the Lagrange multiplier enforcing the volume constraint has a physical interpretation: it represents the pressure difference between solid and liquid, λ * = P = P s − P l , and (2.3) is the Young-Laplace equation relating pressure differences to mean curvature.
Also from δJ = 0 it follows that the following force balance holds along triple lines where three surfaces meet [18]: where ν i is the co-normal to each surface at the triple line and γ i is the corresponding surface energy per unit area. This reduces to the familiar expression for the dihedral angle θ at which two solid-liquid surfaces meet a solid-solid surface (figure 3), It should be noted that the approach taken here, of discretizing the variational problem (minimize J in (2.2)), differs from the approach taken in the classic studies by von Bargen & Waff [5], Cheadle [6] and Nye [7], and also the more recent work by Ghanbarzadeh et al. [10]. The starting point for all of these studies is the statement that the solid-liquid surfaces are of constant mean curvature and meet the solid-solid surfaces at the dihedral angle. These studies are all based on a numerical discretization of the mean curvature of the solid-liquid surfaces, whereas here in discretizing J only areas and volumes are discretized. The approach taken here, working from the variational problem, is essentially identical to the approach pioneered by Beeré [1], although here a more refined discretization is used.
By scaling, the variational problem can be reduced to being a function of just two dimensionless parameters: the porosity φ (the volume fraction of melt) and the dihedral angle θ. For example, if E represents the surface energy contained in one tetrakaidecahedral cell, then a scaled energy can be defined by . All results in this article are presented using scaled variables as functions of φ and θ .

Melt topologies
One of the key features of textural equilibrium is that different kinds of melt topology are possible for different values of the dihedral angle and at different porosities [2,3]. For example, there is the well-known result that, at small porosities, the melt network is only connected along the grain edges for dihedral angles less than 60 • . The problem tackled here allows for a rich variety of melt topologies as a function of porosity and dihedral angle, which are summarized in the regime diagram of figure 4 and discussed below. For some combinations of parameters it is possible to find multiple topologies: each is a local minimum of the energy functional, but not necessarily a global minimum.
The first type of topology, marked 'c' in figure 4, has melt forming a connected network along the grain edges. Examples of such a topology are shown in figure 5, and are broadly similar to those depicted by von Bargen & Waff [5], Cheadle [6] and Nye [7], except that the melt junction does not have tetrahedral symmetry. The second type of topology, 'i' (figure 6), which consists of melt isolated at grain corners, is also broadly similar to previous calculations which assumed tetrahedral symmetry, for which analytical expressions are available [2]. At very large porosities the grain-grain contacts disappear and grains are completely surrounded by melt. This is marked as 'd' for disaggregated in figure 4, where the minimum energy configuration simply consists of spherical grains. 'd' topologies exist for all porosities greater than 1 − π √ 3/8 ≈ 32%. The problem considered here admits additional melt topologies because of the lower degree of symmetry of the quadruple junction. These are depicted in figures 7 and 8. The first of these represents the case 's', where the square faces become wetted but the hexagonal faces do not. The second represents the case 'h', where the hexagonal faces are wetted but the square faces are not.
For low porosities, the dihedral angles for which topologies 'c' and 'i' are found are essentially the same as found for the problems with tetrahedral symmetry: isolated solutions exist only for dihedral angles greater than 60 • . For dihedral angles greater than this, there is a region of overlap between the 'c' and 'i' topologies as porosity increases (the region between the 'pinch-off' and 'wetting' boundaries as described by von Bargen & Waff [5]).
Perhaps the most important new transition in the present study is the transition from 'c' to 's' topologies as porosity increases (i.e. the preferential wetting of the smaller square faces). For low dihedral angles, this occurs close to φ = 0.11, but for larger dihedral angles there is an overlap where both 'c' and 's' topologies exist. The wetting of square faces at around φ = 0.11 is also seen in calculations of wet Kelvin foams [19][20][21][22][23][24]. Wet Kelvin foams are essentially a limiting case of the problem considered here: the case where the dihedral angle approaches 0 • . The network of melt considered here is termed a Plateau border network in the foam literature. The only difference is that a wet Kelvin foam allows for some curvature along what are the grain-grain contacts in the rspa.royalsocietypublishing.org Proc. R. Soc.  In the middle, in pink, is a view of an individual grain. On the right is a view of an individual grain + melt. Corresponding labels give the porosity φ and dihedral angle θ . The first two rows have positive mean curvature; the final row has negative mean curvature.
present model, but such curvature is extremely slight, and unlikely to significantly affect where the transition to wetted square faces occurs.
One curious feature of figure 4 is that for dihedral angles just below 100 • the 's' solutions exist down to very small porosities. If such a solution were realizable, there would be the potential for  percolation of melt at small porosities even for some dihedral angles greater than 60 • . However, as will be seen in the next section, such 's' topologies have higher energy than the 'i' topologies and are thus less likely to be realized.
It should be noted that figure 4 does not provide a map of all possible melt topologies, only a key subset that have been chosen to be investigated. Firstly, the topologies investigated are only those consistent with the symmetries imposed. Topologies with lower degrees of symmetry are possible, e.g. isolated topologies with melt at some quadruple junctions but not others. Moreover, there are additional isolated topologies possible that are consistent with the imposed symmetry but that have not been calculated here. For example, melt can also be isolated in the middle of the grain edges, or at the centres of the grain faces. One can also have some combination of isolated melt in all three places-edges, faces and vertices of the unit cell. However, at the onset of melting it is expected that melt first forms at the quadruple junctions, so it is the case of melt at these junctions that is typically of most interest. Figure 9 plots a scaled surface energy for each of the melt topologies as a function of porosity and dihedral angle. For those regions of φ − θ space that admit multiple topologies, figure 9 can be used to identify which topology has the overall lowest energy (and hence is more likely to be found). For example, the 'h' topologies (which wet the hexagonal faces) always have higher energy than the 's' topologies (which wet the square faces). This is what one might intuitively expect-it is more likely for smaller faces to become wetted than larger faces.

Energy and pressure
The situation is more complicated in the region of overlap between 'c' and 's'. For small dihedral angles, the 'c' topologies have lower energy than the 's' topologies for most of the region of overlap, until close to the boundary where only 's' exists. For larger dihedral angles the 's' topologies are lower energy for most of the range of overlap. Similarly, in the region of overlap between 's' and 'i', for lower dihedral angles the 'i' topologies are lower energy for much of the region of overlap, whereas at higher dihedral angles the 's' topologies are lower energy for the whole region of overlap. Figure 10 plots a scaled version of the pressure difference P between solid and liquid for the various melt topologies. Owing to the Laplace-Young equation (2.3), this is also a plot of a scaled mean curvature. This pressure difference is calculated during the solution of the variational problem as P = λ * is the Lagrange multiplier that enforces the fixed volume constraint. The curves in figure 10 are related to the slopes of the curves for energy in figure 9, because of the relationship which is a consequence of the variational problem. Here, E * represents the energy at equilibrium, and V * l and V * cell the volumes of the liquid and unit cell, respectively. In the scaled (dimensionless) variables used in figures 9 and 10, this relationship becomes Figure 10 shows a number of features that are consistent with previous studies. For example, for the isolated topologies 'i', the curvature is positive for θ < 71 • and negative for θ > 71 • [2]. Interestingly, topologies 'i', 'c' and 's' have zero mean curvature for some combinations of porosity and dihedral angle, and are thus examples of 'minimal surfaces' (surfaces defined as having zero mean curvature).
The pressure difference between the two phases becomes singular as the porosity approaches zero, but the form of singularity varies with the dihedral angle. For isolated topologies the porosity dependence is of the form , plotted as a function of porosity φ for different values of the dihedral angle θ (key on the right). In this plot, and subsequent plots, thin solid lines represent connected topology 'c' , thick solid lines represent square-wetted topology 's' , dotted lines represent isolated topology 'i' and dashed lines represent hex-wetted topology 'h' . The black solid line represents disaggregated topology 'd' (here independent of the dihedral angle). angles, the melt geometry can be closely approximated by tubes along the grain edges, for which P ∝ φ −1/2 (see (B 13) in appendix B).

Effective pressure
The pressure plotted in figure 10 is a measure of the change in energy of the system with liquid content, assuming the cell size remains fixed (and thus as the volume of liquid increases, the  volume of solid must decrease). As pointed out by Park & Yoon [4], in several physical situations what one would like to know is the change in energy with liquid content, but keeping the solid volume (V * s ) fixed. For example, in problems of compaction or sintering one may have a partially molten medium surrounded by a reservoir of liquid, and want to know whether liquid will be drawn in or expelled from the partially molten medium. Park & Yoon [4] quantify this in terms of an 'effective pressure', P e , defined by  [4]. They are scaled such that, as porosity varies, the solid volume remains fixed (rather than the cell size). Scaling is made with respect to the length d 0 , which is defined as the distance between grain centres in the absence of melt with the same volume of solid grain. The area of the unit cell in the absence of melt is given as A cell0 = 3 4 (1 + 2 √ 3)d 2 0 , and the volume as V cell0 = d 3 0 /2. The plot of effective pressure in figure 12 is related to the slope of figure 11 by If the sign of P e is positive, it indicates that there is a driving force sucking liquid into the partially molten material. If negative, the driving force acts to expel liquid. Zero P e represents a stable state (a minimum of the scaled energy plotted in figure 11). The key result of figures 11 and 12 is that topologies for dihedral angles greater than 60 • are always unstable, and the driving force is such that the two phases try to separate (P e is always negative). For dihedral angles less than 60 • there exists a critical porosity at which the effective pressure is zero, and hence the topology is stable.  [4], P e = P e d 0 /γ sl . a porosity less than this critical porosity there will be a tendency of melt to be drawn in, whereas if it is above this critical porosity there will be a tendency for melt to be expelled.
In Park & Yoon's original study, solutions with zero effective pressure could be found for dihedral angles up to 75 • , rather than the 60 • limit found here. One difference between Park & Yoon's study and the present study is that Park & Yoon [4] consider a unit cell taking the shape of a rhombic dodecahedron. However, the discrepancy with the results of this study is likely to have arisen from Park & Yoon's approximation of the grain-melt surfaces by circular arcs, which is not consistent with the surfaces being of constant mean curvature.

Geometrical properties
which measures the relative area of grain-grain contact to total grain area. Figures 14 and 15 show the individual areas of solid-solid and solid-liquid contact, normalized to the area of the unit cell. Contiguity is used as a fundamental variable in a number of micromechanical models, e.g. in determining elastic properties [25] and effective viscosities due to creep [26]. For small porosities, figure 13 shows the same main trends as observed by von Bargen & Waff [5] and Cheadle [6], with contiguity being larger for larger dihedral angles. Contiguity shows a steady decrease with increasing porosity as more of the grain becomes wetted. There is a notable change in slope during the transition from 'c' to 's' topologies, with a shallower slope for the 's' topologies than the 'c' topologies.
In this problem, there are two different kinds of grain-grain contact-those associated with the square faces and those associated with the hexagonal faces. Figure 16 plots the areas of the square faces, and figure 17 that of the hexagonal faces. For the 'c' topologies, the areas of the square faces decrease as porosity increases, but note that they do not vanish at the boundary between the 'c+s' and 's' topologies: thus wetting of the square faces occurs discontinuously from a finite area of square grain-grain contact to a zero area as porosity increases. Only for a  zero dihedral angle does the 'c' topology smoothly go to zero area of the square face as porosity increases. Such discontinuous jumps in behaviour are common in area minimization problems. The classic example of this is in the minimal surface of revolution problem: finding the soap film which minimizes the area between two parallel circular hoops [27]. Beyond a critical separation, the catenoid solution which joins the two hoops breaks down to the Goldschmidt discontinuous solution with two separate films spanning the two circles. For small dihedral angles, there is a close agreement between these results and those for wet Kelvin foams: hexagonal and square areas reported in figure 5 figure 18). For small porosities, and small dihedral angles, figure 18 shows that this value is approached, reflecting the fact that for such porosities and dihedral angles the melt geometry is well approximated by tubes. At small porosities, the mid-edge channel area decreases with increasing dihedral angle, reflecting the fact that more melt resides near the grain vertices than resides near the grain edges. The opposite trend is seen for the 's' topologies at large porosity, where the larger dihedral angles have larger channel areas.

10°2
0°3 igure 18. Scaled cross-sectional channel area (sectioned at the mid-grain edge), scaled as A ch /(d 2 φ). Isolated topologies 'i' have no melt at mid-grain edges and thus have zero area and are not shown. The black line labelled 'u' on the y-axis shows the expected channel area for melt tubes of uniform cross section along the grain edges.

Permeability
The effective permeabilities of the melt topologies are shown in figure 19, and details of their calculation can be found in appendix C. Permeabilities are scaled in figure 19 by d 2 φ 2 , based on the expected behaviour for low porosities (see below). Most of the main features of figure 19 can be understood in the context of simpler models. For small porosities, and small dihedral angles, the melt geometry for topology 'c' takes the form of tubes along the grain edges. A very simple model of permeability was derived by Frank [  section lying on the edges of a tetrakaidecahedron. For this model, the permeability is cross section. For a cross-sectional shape more appropriate for low dihedral angles (three circular arcs meeting at zero dihedral angle), the melt flux for a given cross-sectional area and pressure gradient is a factor of 0.505 lower. A modified version of Frank's formula is then and this formula is marked as 'u' on the y-axis of figure 19. The equation above explains the limiting behaviour seen in figure 19 for low porosities and low dihedral angles, which tend to finite values of k/(φ 2 d 2 ) as φ → 0, and whose geometries can be well described by tubes along the grain edges with approximately uniform cross section. For larger dihedral angles the cross sections become much less uniform, with more melt residing at the vertices of grains rather than the edges. The effect of this is to reduce the permeability. von Bargen & Waff [5] provide an approximate formula, which is in good agreement with the results here for small porosities and dihedral angles around 40 • . Such a formula gives the permeability to within 15% for all porosities less than 2%. The expression quoted by Cheadle [6] (k = φ 2 d 2 /3000) seems to underestimate the permeability, and the permeabilities calculated here are much closer to the results of von Bargen & Waff [5]. For larger porosities (greater than 10%), lines on figure 19 are roughly linear, which reflects the fact that at large porosities permeability scales as k ∝ d 2 φ 3 . The φ 3 behaviour can also be understood from simpler models, and is the expected behaviour for melt lying as thin sheets along the grain boundaries. A reasonable approximation for 0.1 < φ < 0.3 and the 's' topologies is which lies within about 25% of the calculated values for dihedral angles θ < 60 • . Indeed, it should be noted that the relative range in permeability is very small when θ varies for porosities above 10%: a factor of 2 at most. It is interesting to note that the behaviour of permeability with dihedral angle differs at higher porosities: for high porosity, the larger the dihedral angle, the greater is the permeability; although the magnitude of this effect is rather modest. Permeability data are often parameterized with a single law of the form k = φ n d 2 /C (e.g. [30,31]). Note that such a parametrization in terms of a single C and n value would not be able to describe all the subtleties seen in figure 19, and particularly the shift from k ∝ d 2 φ 2 at low porosities to k ∝ d 2 φ 3 at larger porosities.

Discussion
As remarked in the Introduction, the textural equilibrium geometries calculated here have also been calculated by Ghanbarzadeh et al. [10][11][12][13][14], although these authors consider a more limited range of porosities and dihedral angles. Ghanbarzadeh et al. [11] introduce a novel approach to calculating textural equilibrium geometries based on the level-set method. In this approach, the interface between the solid and liquid is represented by the level set of a function, and the interface evolves at a velocity proportional to the surface Laplacian of mean curvature. Over time, the surface should approach a state of constant mean curvature. A partial differential equation determines the evolution of the level-set function, and additional terms are added to the PDE to enforce the dihedral angle constraint. Ghanbarzadeh et al. [11] discretize the PDE using high-order finite differences.
One of the advantages of Ghanbarzadeh et al.'s approach is that it allows a great deal of geometric flexibility, including the ability to solve for many grains at once with different grain shapes. As such, they do not explicitly impose symmetry on the solution (as is done here). Their method is considerably more computationally expensive than the approach taken here (which can solve for the melt geometries in seconds). In broad terms, many of the topologies calculated by Ghanbarzadeh et al. [11] are similar to those calculated here, but there are some notable discrepancies. One of the most important discrepancies concerns the wetting of the square faces. Ghanbarzadeh et al. [10] state that 'In contrast with prevailing assumptions, the smaller square grain boundaries become fully wetted at θ = 10 • and φ ≥ 5%.' As can be seen in figure 4, this study places the transition to wetted square faces at closer to φ = 11% for θ = 10 • : a factor of 2 difference. As remarked earlier, transitioning at φ = 11% is consistent with results for wet Kelvin foams [19,22,23].
Another discrepancy is in the calculated mean curvatures of the melt topologies. Hesse et al. [32] presented a plot of mean curvature against porosity for isolated 'i' and connected 'c' topologies for a dihedral angle θ = 90 • . In theory, this plot should be similar to the curvature information shown in figure 10, but in fact bears little resemblance. First, Hesse et al. [32] showed 'c' topologies for porosities as low as 1.5%, whereas the suggested lower bound for 'c' topologies here is around 4% (very similar to the pinch-off boundary in fig. 7 of von Bargen & Waff [5]). Moreover, the behaviour of the mean curvature with porosity for the isolated 'i' topologies shown by Hesse et al. [32] is very different from that seen here. The shape of the isolated topologies is independent of porosity, apart from a scaling of the coordinates [2]. Thus mean curvature H should be proportional to φ −1/3 for the 'i' topologies, and be singular at zero porosity (as indeed seen in figure 10). The curve shown by Hesse et al. [32] does not show this behaviour, and instead appears to curve the wrong way, and seems to approach a finite value for small porosity.
These discrepancies matter, because they lead to different predictions about upscaled physical properties like permeability. Indeed, the permeabilities that have been reported by Ghanbarzadeh et al. [13,14] differ from those calculated here by a significant margin, in some cases by almost an order of magnitude. For example, for a dihedral angle of 10 • , Ghanbarzadeh et al. [13,14] report that their permeability data can be fitted well by an expression of the form k = d 2 φ 2.6 /595.66. For a porosity of φ = 0.01 this implies k/(d 2 φ 2 ) ≈ 1 × 10 −4 , whereas here a value of k/(d 2 φ 2 ) ≈ 7 × 10 −4 has been estimated ( figure 19). Thus, Ghanbarzadeh et al. [13,14] seem to underestimate the permeability by a factor of 7 for these parameter values. For a porosity of φ = 0.1, Ghanbarzadeh et al.'s expression implies k/(d 2 φ 2 ) ≈ 4 × 10 −4 , whereas here k/(d 2 φ 2 ) ≈ 1.4 × 10 −3 has been estimated, more than a factor of 3 greater. Perhaps the most probable reason for the discrepancies is that the simulations of Ghanbarzadeh et al. [11] are under-resolved, and struggle to accurately capture small-scale variations in the geometry, e.g. with cusps at a small dihedral angle or small radii of curvature at low porosities. The accuracy of the solutions obtained in this study is discussed in appendix A.
Another point of difference between this study and that of Ghanbarzadeh et al. is that Ghanbarzadeh et al.'s study essentially aims to produce surfaces of constant mean curvature compatible with the dihedral angle constraints. Constant mean curvature is a necessary, but not a sufficient, condition for minimum energy. Constant mean curvature guarantees that the energy is extremized, but does not say whether the surface is a minimum, maximum or a saddle point of the energy. Indeed, it is well known that some solutions of the Euler-Lagrange equations for area minimization problems can be unstable, as in the minimal surface area of revolution problem mentioned earlier. Hence it is possible that the scheme proposed by Ghanbarzadeh et al. produces melt topologies that are not minimum energy, and hence would not be found by the approach taken here.
The existence of multiple melt topologies has important consequences for upscaled physical properties (e.g. permeability, electrical conductivity, effective viscosities). First, transitions in such properties can be discontinuous on varying parameters, as the underlying transitions in topology are discontinuous. Moreover, as remarked by von Bargen & Waff [5], and more recently by Hesse and co-workers [14,32], there is the potential for hysteresis in topology and hence hysteresis in upscaled physical properties. For example, suppose the dihedral angle is 40 • and porosity is increased slowly from zero such that textural equilibrium is maintained ( figure 4). The topology would remain as type 'c' until a porosity of 12% when a discontinuous transition to the type 's' topology would occur (wetting of the square faces), when the 'c' state is no longer stable. If porosity were then slowly reduced, the topology would stay as type 's' until a porosity of 9% when a discontinuous transition back to 'c' would occur, as the grain-melt surfaces begin to touch on the square faces (de-wetting of the square faces). Similar hysteresis can occur in other parts of the regime diagram, for example the 'i' to 'c' transition discussed by von Bargen & Waff [5].
There are several natural ways in which the present work can be extended. In the same way that this study presents a model with a lower degree of symmetry than the studies by von Bargen & Waff [5], Cheadle [6] and Nye [7], it would be fruitful to examine the consequences of relaxing more of the symmetry assumptions. Symmetry plays a crucial role in constraining the solutions produced here, and topologies which are stable under the symmetry constraints imposed here may not be stable if some of these symmetry constraints are relaxed.
In this study, a tetrakaidecahedral unit cell was chosen, but it would be straightforward to repeat the same calculations for other choices of unit cell, such as the rhombic dodecahedral unit cell used by Park & Yoon [4]. In a tessellation of tetrakaidecahedrons, at each vertex four edges meet, and each vertex is thus said to have a coordination number of 4. In a tessellation of rhombic dodecahedrons, half of the vertices have a coordination number of 4, and half have a coordination number of 6. X-ray tomographic imaging suggests that the predominant coordination number in olivine-basalt systems is 4 [33], so the tetrakaidecahedral unit cell is the more appropriate geometry for the partially molten rocks of the Earth's mantle.
To simplify the calculations, grain boundaries were forced to be planar, and forced to lie on certain symmetry planes. Grain boundaries are in general not perfectly planar, and these constraints could be relaxed in future work. Indeed it is the curvature of grain boundaries that allows for grain growth and coarsening, an effect that is suppressed in the present calculations by the assumptions that have been made. Another natural extension is to consider the anisotropy in surface energy, both between grains and between the grains and the melt. Anisotropy can lead to faceted solid-liquid and solid-solid interfaces (e.g. [34]) that will present challenges for numerical simulation. Another challenge for the future is to move beyond textural equilibrium, to study the interplay between deformation and melt network geometry.

Conclusion
The aim of this study has been to produce a simple reference model which describes the geometry of melt within a partially molten material. This model, based around a tessellation of tetrakaidecahedral unit cells, has an advantage over most previous models in that the melt geometry is compatible with a space-filling tessellation of grains. A lot of the fundamental behaviour of the melt topologies calculated here agrees with the classic models of von Bargen & Waff [5], Cheadle [6] and Nye [7], particularly at small porosities. The main differences occur at higher porosities, where topologies that wet the square faces of the tetrakaidecahedral grains exist.
These geometries form a starting point for the calculation of upscaled physical properties which depend on the geometry of melt at the grain scale. One example of such a property has been presented here, namely the permeability, but the intent is to describe a wider range of physical properties in future articles. In particular, calculations on the effective creep properties will be presented, building on the simpler grain-scale models of Takei

(c) Different topologies
The different topologies explored in this contribution have different problem set-ups, essentially due to the removal of one or more of the bounding planes from the connected 'c' topology problem. For example, to find 's' solutions that wet the square faces, the bounding plane corresponding to the square triple line was removed. To find isolated 'i' solutions, the two planes associated with truncation were removed. To find hex-wetted 'h' solutions, the plane associated with the hexagon triple line was removed.

(d) Accuracy of solutions
The accuracy of the numerical solutions has been checked by comparing them against analytical solutions for combinations of parameters for which such solutions exist (see appendix B). For example, for an 's' type solution with a dihedral angle of 60 • that consists of spherical surfaces, the numerical mean curvature differs from the analytical expression (B 8) by only 0.004%. As the spacing between mesh edges has been kept fixed, the solutions at low porosity are liable to be less accurate than those at higher porosity as fewer mesh nodes are used. However, the good comparison with the expected behaviour at low porosity suggests they are well resolved. between different topologies were determined by two different methods depending on the nature of the boundary. The first type of boundary occurs when the topology assumed in the calculation is inconsistent with the tessellation. For example, when calculating the isolated topologies, there is a point at which the individual isolated pockets of melt meet as porosity increases. This can be determined by finding the porosity at which the extent of the isolated pocket is such that it extends to the middle of a grain edge. Similarly, the lower boundary of the 's' topology can be determined by when the grain-melt surface near the square face first touches the boundary of the unit cell. This type of boundary is straightforward to determine accurately, and can be obtained through a simple interpolation of results at different porosities. The second type of boundary is more difficult to constrain. This occurs when the topology in question becomes unstable as the parameters change. An example of this type of boundary is where 'c' topologies no longer exist. This type of boundary can be identified by the lowest eigenvalue of the Hessian matrix [16], as exemplified by figure 20. As porosity increases, the lowest eigenvalue decreases until at some point it reaches zero and the surface is no longer a minimum of the energy, but a saddle point. Following [16], the boundaries of such topologies have been estimated based on linearly extrapolating to the point at which the lowest eigenvalue reaches zero. As the method is based on extrapolation, it is inherently somewhat inaccurate. Additional simulations have been run in the vicinity of boundaries to improve the accuracy (see additional points near 12% in figure 20). Such boundaries are expected to be accurate to within at least ±0.5% and ±5 • .

Appendix B. Analytical and approximate solutions (a) Spherical surfaces
In certain special cases, the melt geometry can be calculated analytically, and these cases form useful tests of the numerics. These analytical cases arise when the solid-liquid interface is spherical. The first of these cases is for the disaggregated topologies 'd', where the grains are isolated spheres. The energy and curvature are then given by E = γ sl d 2 π 1/3 (3(1 − φ)) 2/3 (B 1) and The second case arises for the 'i' topologies, when the dihedral angle is θ = 180 • and the melt forms spherical bubbles at the grain corners. Quantities of interest in this case are theory (e.g. [28,35]). In dimensionless form, the cell problem to solve is ∇ 2 u j = ∇P j − e j (C 1) and ∇ · u j = 0 ( C 2 ) for velocities u j and pressures P j ; e j is the jth unit vector (j = 1, 2, 3). No-slip (u j = 0) boundary conditions are applied to the solid-liquid boundaries, and appropriate periodic boundary conditions are applied on the boundaries of the unit cell. The permeability tensor is then given by where u ji is the ith component of velocity u j and V is the volume of the unit cell. Owing to the cubic symmetry, the permeability tensor must be isotropic, and hence it suffices to solve the Stokes flow problem only for j = 1. A tetrahedral mesh of the liquid domain for the unit cell was generated with mesh edge size equal to that used for the triangulation of the solid-liquid surface. Equations (C 1) and (C 2) were solved numerically using Taylor-Hood finite elements (second-order velocity, first-order pressure) with the FENICS software [36,37].