## Abstract

A model for the formation of granitoid systems is developed involving melt production spatially below a rising isotherm that defines melt initiation. Production of the melt volumes necessary to form granitoid complexes within 10^{4}–10^{7} years demands control of the isotherm velocity by melt advection. This velocity is one control on the melt flux generated spatially just above the melt isotherm, which is the control valve for the behaviour of the complete granitoid system. Melt transport occurs in conduits initiated as sheets or tubes comprising melt inclusions arising from Gurson–Tvergaard constitutive behaviour. Such conduits appear as leucosomes parallel to lineations and foliations, and ductile and brittle dykes. The melt flux generated at the melt isotherm controls the position of the melt solidus isotherm and hence the physical height of the Transport/Emplacement Zone. A conduit width-selection process, driven by changes in melt viscosity and constitutive behaviour, operates within the Transport Zone to progressively increase the width of apertures upwards. Melt can also be driven horizontally by gradients in topography; these horizontal fluxes can be similar in magnitude to vertical fluxes. Fluxes induced by deformation can compete with both buoyancy and topographic-driven flow over all length scales and results locally in transient ‘ponds’ of melt. Pluton emplacement is controlled by the transition in constitutive behaviour of the melt/magma from elastic–viscous at high temperatures to elastic–plastic–viscous approaching the melt solidus enabling finite thickness plutons to develop. The system involves coupled feedback processes that grow at the expense of heat supplied to the system and compete with melt advection. The result is that limits are placed on the size and time scale of the system. Optimal characteristics of the system coincide with a state of maximum entropy production rate.

## 1. Introduction

The Earth cools by transferring heat from its interior to the base of its outer thermal boundary layer (the lithosphere) by convective mechanisms. Transfer of heat through the lithosphere is by three mechanisms, namely (i) conduction in the solid state, (ii) advection of heat upwards through the migration of melts, and (iii) advection of heat upwards through the migration of both hydrous and non-hydrous (for example CO_{2}-, CH_{4}-bearing) fluids. Ultimately the heat is radiated into space through the atmosphere. Conduction of heat is perhaps the dominant mode of heat transfer through the lithosphere. This is a steady-state process on time scales, *τ*, measured by *τ*=*L*^{2}/*κ*_{thermal} where *L* is the thickness of the crust or lithosphere and *κ*_{thermal} is the thermal diffusivity (taken here to be 10^{−6} m^{2} s^{−1}). Thus *τ* varies from 50 Myr for a crust 40 km thick to approximately 260 Myr for a lithosphere 120 km thick. These time scales are long compared with those (approx. 10^{4}–10^{7} years) proposed for the development of granitoid complexes (Brown 2001, 2010; Coleman *et al*. 2004; Glazner *et al*. 2004). This compares with longer time scales (at least 30 Myr) for the duration of high-grade metamorphism (Rubatto *et al*. 2001; Willigers *et al*. 2001). Thus it seems that processes other than steady-state heat conduction are responsible for the formation of granitoid systems and these processes are short lived compared with the lifetime of an orogen.

Advection of heat by the upward migration of melts and/or fluids is a more efficient means of transferring heat out of the Earth than conduction. The upward migration of melts can lead to the formation of large igneous complexes whereas the upward migration of fluids leads to the formation of giant hydrothermal ore deposits under special conditions. Both processes transfer heat upwards faster than conduction. In some instances the two processes are closely coupled but it is the intent in this paper to consider only the migration of melts and then only those associated with the development of granitoid systems such as occur in the Himalayan orogenic zone (Hodges 1998) or in the Tasman Orogen in Eastern Australia (Vernon & Clarke 2008).

All of the processes mentioned above represent non-equilibrium states and so models of granitoid systems need to be developed using the principles of modern continuum thermodynamics (Zeigler 1983*a*; Coussy 2004; Houlsby & Puzrin 2006). This paper is a beginning on a path towards such models. One of the tenets of modern non-equilibrium thermodynamics is that systems not constrained to be at steady state will evolve to that state which tends to maximize entropy production. Such a principle has proved invaluable in explaining aspects of the Earth’s climate (Paltridge 1975, 1978, 2001), the oceanic general circulation (Shimokawa & Ozawa 2002), and the fundamental mechanical, hydraulic and chemical behaviour of solids (Ziegler 1983*a*,*b*; Coussy 1995, 2004; Houlsby & Puzrin 2006). The principle has a basis in statistical mechanics (Dewar 2005) where it has been shown that a state of maximum entropy production rate is the most probable state for large systems not at equilibrium and not constrained to be at steady state.

In this paper, the aim is to show that granitoid systems in the crust of the Earth evolve to a state that maximizes the entropy production rate through a complex interaction between the production of melts, stable and unstable flow and retrograde reactions within the Melt Zone. It is well documented (Prigogine 1955; Kondepudi & Prigogine 1998) that thermodynamically linear steady-state processes such as thermal conduction minimize the entropy production rate. The principle of maximum entropy production rate sets limits on the rates (and hence the size) at which systems can evolve once unsteady processes begin to operate.

## 2. Entropy production rate in granitoid systems

The crust of the Earth is an open dissipative system that undergoes elastic and inelastic deformation while it receives thermal and mass fluxes from the mantle of the Earth. Many granitoid systems originate from fluid-absent melting of crustal (mainly meta-sedimentary) materials with no addition of mass from the underlying lithosphere (Brown 2001). Other models (Collins *et al*. 2000) involve addition of mafic melts from the mantle and/or addition of hydrous fluids either from a dehydrating subducting slab or some other source in the mantle or from inside the crust (Brown 2001; Berger *et al*. 2008). The crust also loses heat and mass to the hydrosphere. The mass fluxes consist of both volcanics that leave the crust and form deposits at its interface with the hydrosphere, and fluids such as H_{2}O, CO_{2} and SO_{2} that are exsolved from the melt during crystallization. While these processes are operating, the system is not at equilibrium and for some of its life the system is evolving rapidly and is not at steady state. In such systems it appears that the rate of entropy production is maximized (the principle of maximum entropy production rate; Ziegler (1983*a*)) and so it is therefore of interest to explore the implications of this principle for the evolution of the crust.

In this paper, we consider only the formation and migration of granitic melts and their emplacement as plutons in an attempt to understand more about the systematics of magma accumulation in the crust. In a subsequent paper we will consider the addition of other fluids in an attempt to understand more about the systematics of giant hydrothermal ore body formation. The lithosphere can lose mass to the hydrosphere through igneous eruptions; again we neglect this issue here. We also neglect the addition of fluids to the lithosphere from the hydrosphere.

The dissipative processes that operate in this system and that therefore contribute to or detract from entropy production are (i) inelastic deformation, (ii) mineral reactions, both those that produce melt and retrograde reactions that produce hydrous phases, (iii) thermal transport by conduction and by advection with the moving magma, (iv) the flux of melt across gradients in temperatures within the Melt Zone, (v) introduction of melt into the stressed crust; this is akin to adding a chemical component into a stressed solid and so in principle is conceptually the same as defining a ‘chemical potential’ for the melt, (vi) the transport of mass and/or heat across the top and bottom of the crust, and (vii) the internal generation of heat due to radioactive decay of U, Th and K. In order for the rate of entropy production to reach a maximum, and not increase indefinitely, some of the processes mentioned above must evolve to compete with others. Some of these processes such as inelastic deformation, non-hydrating silicate reactions, melting and introduction of melt into a stressed crust are endothermic. Others such as hydrating retrograde reactions and crystallization of melts are exothermic. This means that in the granitoid system as a whole various processes compete for the supply and consumption of heat (that is, entropy). The essential aim of this paper is to examine the contributions these processes make to entropy production rates in the system and how this influences the spatial distribution and total volume of granitoid intrusions in the crust.

The specific questions we seek to address in this paper are as follows: (i) What controls the total volume of crustally derived magma emplaced into the middle-to-upper crust? (ii) What controls the level at which this magma accumulates in the crust? (iii) Is there an optimal set of conditions that produces giant granitoid complexes such as occur in the Ashuanipi sub province, Quebec (Gurnina & Sawyer 2003) and the Lachlan Fold Belt of Eastern Australia (Vernon & Clarke 2008)? (iv) What sets the upper limit to the production and size of such systems? The fundamental hypothesis is that all of these questions are answered if one can understand the conditions under which the system evolves to maximize the rate of entropy production. Before proceeding to examine that hypothesis we need to define the evolution of granitoid systems in a way that emphasizes the broad processes involved; we remain at this macro-level throughout and delve into the detailed mechanisms and processes that operate only where necessary. The arguments developed here are meant to be quite general. Application to specific examples requires knowledge of local parameters, in particular the thermal history together with the chemical and mineralogical composition of rocks in the anatectic zone. Even so, it is hoped that the general results presented here will encourage field geologists to examine granitoid systems with a new set of eyes.

The system is illustrated in figure 1 (table 1) and comprises the following. (i) The mantle of the Earth, which is assumed to be at steady state and losing heat into the crust of the Earth, via the heat flux, *q*_{IN}. (ii) The crust of the Earth, which is being deformed via a stress *σ*_{ij}. The crust is not at steady state and is dissipating energy arising from deformation, melting, and transport and emplacement of melt. All of this energy is expended within the crust. The heat flux from the top of the crust is *q*_{OUT}. There is no mass flux out of the crust into the atmosphere. (iii) The atmosphere, which is assumed to be at steady state with a temperature at its base of *T*_{0}.

The crust comprises the following. (i) The Melt Zone where latent heat is consumed. Its top is defined by a moving isotherm that defines the beginning of melting; this is the anatectic front. (ii) The Melt Transport Zone. Its top is defined by a temperature such that the constitutive properties of the magma change from elastic–viscous to elastic–plastic–viscous. This has two subzones as discussed in the text. (iii) The Melt Emplacement Zone. Latent heat is completely released here and its top is defined by the melt solidus isotherm; also, at this level, the melt pressure gradient drops below lithostatic. (iv) The Cap Zone, whose top is defined by the surface of the Earth. Radiogenic heat sources exist in each zone.

## 3. A model for granitoid systems: an overview

We present below an overview of the model for granitoid systems developed in detail in the remainder of this paper. The term ‘granitoid system’ is meant to refer to the body of rock starting at the base of the crust and comprising (i) a melt source region, (ii) a melt transport region, and (iii) a melt emplacement region. The melting process begins as rising isotherms first induce a temperature corresponding to the wet solidus (if any initial H_{2}O is present) and then progressively to each of the generalized H_{2}O-absent melting reactions muscovite-breakdown melting, biotite-breakdown melting and hornblende-breakdown melting. Examples of specific H_{2}O-absent reactions include the following:

Reaction I: Ms+Ab+Qtz→Kfs+Als+L (Peto 1976);

Reaction II: Ms+Pl+Qtz→Kfs+Sil+Bt+L (Patino Douce & Harris 1998);

Reaction III: Bt + Pl + Als + Qtz → Grt + Kfs + L (LeBreton & Thompson 1988);

Reaction IV: Bt+Pl+Qtz→Opx+Grt+Kfs+L (Vielzeuf & Montel 1994).

In most instances, in the absence of H_{2}O, muscovite breakdown melting reactions will be encountered first, but we will show that the other reactions involving biotite and hornblende are important in controlling the thermal structure of the Melt Zone and in controlling the supply of melt to the emplacement region. Whether they are encountered or not depends on the *P*–*T*–*t* paths characterizing the metamorphic history of the system and on the chemical composition of the rocks involved. We show that the melt pressure at the top of the anatectic front depends on the volume concentration of H_{2}O and/or the hydrate phase, the rate at which the melting isotherm moves upwards and the volume of melt produced by a particular reaction at that isotherm. This isotherm is the fundamental flux and pressure control valve in the system that influences the entire history of melt transport and emplacement higher in the system. Once melt is generated it must move upwards driven by buoyancy, although local effects arising from deformation may inhibit or enhance this movement. In addition, topography at the surface of the Earth can drive significant horizontal flow. The flux generated at the moving isotherm is relentless and controls the permeability distribution and overall behaviour of the complete system. The isotherm representing first melting rises slowly if conduction is the only mechanism operating so that vertical advection of heat by rising melt in the Melt Zone is necessary to produce granitoid systems rapidly (less than 10^{7} years).

Here we propose a process where melting is thermodynamically coupled to viscous solid-state deformation in the Melt Zone such that both strain- and strain-rate softening result for the mass as a whole. This produces the meso-scale structures that one typically observes in migmatite terrains. These concepts are further developed in Hobbs *et al*. (in press) where the coupling between deformation and mineral reactions, including melting, is considered in detail. The constitutive model for that region of the crust where the constitutive behaviour is strongly rate dependent is modelled as a Gurson–Tvergaard material (Gurson 1977; Tvergaard 1989; Regenauer-Lieb 1999; Nahshon & Hutchinson 2008). This means that although the material can behave as a rate-dependent ‘power-law’ material it can also accommodate the formation, growth and coalescence of melt inclusions. Given the appropriate stress conditions and initial melt inclusion population (size, age and spatial distributions), the melt inclusions can coalesce as sheets (both parallel and inclined to initial foliations) that enable melt transport in the form of ‘ductile dykes’. Within the Melt Zone these processes are self-enhancing and ultimately result in vertical advection of heat that accelerates the upward velocity of the melting isotherm and hence influences the melt flux-pressure control valve at the top of the Melt Zone. The initiation of heat advection by upward-moving melt in the melting regime is a fundamental step in accelerating the upward migration of isotherms and hence in allowing the rapid production of granitoid complexes. The subsequent initiation of reactions involving biotite-breakdown melting and, possibly, hornblende-breakdown melting as the system becomes hotter enhances these processes even further with different reactions operating at different levels in the Melt Zone.

The development of melt inclusion sheets and ductile cracks in the transport zone above the reaction isotherm enables the melt to propagate upwards at a rate that is governed by the melt flux developed at the top of the Melt Zone. There is a critical upward velocity of melt that enables heat to be advected fast enough to allow the isotherm corresponding to the melt solidus to be raised to high levels in the crust on the time scales of 10^{4}–10^{7} years so the melt can be transported without solidifying.

Two parts to this transport zone can exist: a lower one where the melt is still in local thermal equilibrium with the surrounding rocks and transport is in closely spaced ductile fractures, and an upper one where the melt is hotter than the surrounding rocks and transport is in widely spaced fractures. A conduit selection process operates in this zone to grow progressively wider conduits as the melt moves upwards. The lowest of these subzones (spatially) is at a temperature corresponding to supra-solidus conditions for the melt but sub-solidus for the host rocks whereas the uppermost subzone is at temperatures corresponding to sub-solidus conditions for the host rocks and is close to the melt solidus but the melt is losing latent heat. The transition from one part of the transport zone to another is controlled by a decrease in magma viscosity as it cools; this viscosity evolution drives the conduit size-selection process and ultimately results in wide conduits that are the feeder dykes for plutons.

Variations in melt pressure at the anatectic front diffuse through the system on a time scale measured by *L*^{2}/*κ*_{melt} where *L* is a length scale and *κ*_{melt} is the pressure diffusivity of the melt whereas temperature relaxation resulting from a pressure decrease at the anatectic front diffuses through the system on a time scale measured by *L*^{2}/*κ*_{thermal}. Since *κ*_{melt}≈4*κ*_{thermal} pressure fluctuations can occur in the system with little short-term influence on the temperature. Ultimately both the melt pressure gradient and the temperature gradient begin to decrease due to either a decrease in heat supply or an absence of suitable rock compositions in the melting region. If the system is tuned such that the pressure gradient high in the system decreases below the lithostatic gradient before, or as, the melt temperature drops below the solidus, then the melt can accumulate in structurally controlled sites that can still open driven by the melt flux. This results in a rapid decrease in temperature and the melt finally solidifies. The final emplacement site is a trade-off between the waning melt pressure and temperature gradients. However, melt accumulation in the form of plutons cannot occur unless the constitutive behaviour of the cooling magma changes from elastic–viscous (typical of fluids) to elastic–plastic–viscous (typical of solids). The progressive initiation of muscovite-, biotite- and hornblende-breakdown reactions in rocks with different fertility and, therefore, melt production rates at any given *T* for fixed *P* means that different pulses of melt can be added to the system resulting from different melt fluxes at the top of the Melt Zone. This leads to the development of stacked or composite plutons (Paterson & Vernon 1995; Glazner *et al*. 2004). There is a critical value of the melt flux (for a given thickness of Melt Zone) at the top of the Melt Zone that results in a maximum for the rate of entropy production. If this flux increases above this critical value, instabilities such as convection within the Melt Zone initiate that lead to smaller entropy production rates. The dramatic increase in melt viscosity associated with the resulting retrograde reactions at the anatectic front is a fundamental step in reducing the rate of entropy production. This critical melt flux at the anatectic front sets an upper limit to the size of the Melt Zone and, hence, of the complete granitoid complex. These conditions represent a state of maximum entropy production rate.

## 4. The mechanics of the heat source region

Melting at the base of the crust arises from four sources that may or may not be coupled. These are the following. (i) An increase in heat flux arising from some thermal event in the lithosphere such as delamination, underplating or processes arising from slab dehydration. We are not concerned with the details of such events; we require only that these events result in an increase in thermal flux through the base of the crust. (ii) The introduction of hot, mafic igneous material from the mantle. Although this undoubtedly leads to complications in melt viscosities such as mixing or mingling of mafic and granitic melts occurs we will neglect such complications and treat the introduction of such mafic material as essentially a heat input into the crust from the mantle. (iii) The introduction of hot fluids such as H_{2}O and CO_{2}. The processes for achieving this are the same as will be discussed later for the migration of melt through the lower part of the Transport Zone. The outcome is advection of heat with the moving fluid and presumably some form of metasomatism of the rocks or incorporation of H_{2}O in the lower crust prior to or during melting. (iv) An increase in the thickness of a crust containing radiogenic heat production elements, in particular U, Th and K arising from overthrusting, folding or sedimentation (England & Thompson 1984). The ultimate outcome of these processes is an increase in the upward speed of isotherms in the crust.

## 5. The mechanics of the melt zone

The top of the Melt Zone is defined by an isotherm that propagates upwards from the base of the crust and corresponds to the initiation of anatexis within the crust governed by the local chemical composition (including H_{2}O content) and pressure. During orogenesis the Melt Zone is also deforming and topography is developing at the surface of the Earth and so melt pressure gradients arise from a combination of the following. (i) Buoyancy derived from contrasts in density between melt and the host rocks. This results in a component of upward flow. (ii) Inelastic deformation that produces gradients in mean stress. These gradients exist on the scale of the structures formed during deformation but the increase in mean stress arising from deformation can never exceed a small proportion (say 10–20%) of the lithostatic pressure without inducing ductile or brittle fracture. This limits the magnitude of the increase in mean stress (at say 20 km depth) to about 50 MPa. Thus melt may be driven from high mean stress regions in the inner hinges of folds towards low mean stress regions such as boudin necks provided these mean stress gradients coincide with melt hydraulic potential gradients. (iii) Variations in lithostatic pressure in the Melt Zone arising from variations in topography at the surface of the Earth. These variations produce horizontal gradients in melt pressure inducing horizontal flow.

The relative contributions of these three processes are indicated in table 2. The buoyancy forces are independent of scale and are fixed by the difference in density between the melt and the host rock. The hydraulic head gradients arising from deformation depend on the scale and the orientation of the pressure gradient. At the scale of 1 km or less such gradients can readily compete with or enhance the buoyancy gradient and even at larger length scales are likely to be important. If the deformation-induced gradients compete with the buoyancy-induced gradients then ponding of magma can occur on the scale (vertically) of 1–5 km, which depends solely on the scale of structures developed in the Melt Zone. Specific examples of such ponding processes are discussed by Connolly & Podladchikov (2004) and Hobbs *et al*. (2004). The ponding effect is not due to the existence of a permeability barrier and arises solely from the hydraulic potential gradient becoming zero over some height due to the deformation-induced gradient offsetting the buoyancy-induced gradient. As such the ponding is transient since it relies on continued deformation and constant temperature and must migrate as the deformation and temperature change. Note from table 2 that horizontal flow arising from topographically driven flow can be close in magnitude to the vertical flow induced by buoyancy or can exceed vertical flow if the horizontal permeability is high enough. This is consistent with the melt flow scenarios proposed by Brown & Solar (1998*b*) for a granite complex in NE, USA, and Weinberg & Searle (1999) and Leitch & Weinberg (2002) for the Himalayan leucogranites.

The upward migration of melt is governed by the following. (i) The speed of the upward migration of the isotherm corresponding to initial melting. (ii) The volume rate of melt production from the relevant melt-reaction. (iii) The fertility or volume concentration of components that can contribute to melting in the rock type at the anatectic front. These statements assume that only one melt-reaction is initiated in the Melt Zone and that melt-reaction is defined only at a particular temperature and pressure but can be generalized to systems where melting occurs over a range of *P* and *T*. This simple situation is not what actually occurs, and, as the isotherms rise, different reactions such as those labelled I–IV are progressively initiated. Moreover, melting for a particular chemical composition does not occur at a specific temperature for a prescribed reaction but takes place over a range of *P*–*T* conditions. An example is described by Vielzeuf & Montel (1994). Thus the Melt Zone is a dynamic region where more than one anatectic front exists within the zone and these fronts move upwards as the system continues to evolve. For simplicity, throughout the remainder of this paper we will refer to the top of the Melt Zone as the anatectic front (figure 1).

In what follows, we first assume that only one reaction takes place and that melting occurs at a specific *P* and *T* with no subsequent melting over a range of *P* and *T*. We then proceed to generalize this simple model.

### (a) Production of melt and the pressure control valve

Following Phillips (1991) we consider a reaction
where *S*_{1} is an assemblage bearing an initial volume concentration, *s*_{0}, of a hydrous phase such as muscovite, biotite or amphibole and *S*_{2} is an anhydrous assemblage. In this reaction, *s* cubic metres of the hydrous phase per cubic metre of the assemblage *S*_{1} decrease with time as they are replaced by *S*_{2}
If *T*_{e} is the equilibrium temperature at which the reaction takes place, for *T*<*T*_{e}, *Q*_{s}=0. When *T*≥*T*_{e} the rate of disappearance of *S*_{1} is proportional to *s* and, near equilibrium, to (*T*−*T*_{e}), so that *Q*_{s} can be written as
where *γ* is the reaction rate. Thus,
and
If the melting isotherm is progressing upward through the rock mass with a speed *U*, then, at any point near the melting isotherm,
where *t* is the time since the melting isotherm passed through that point and grad *T* (more than 0) is the negative temperature gradient in the direction of *U*. Thus,
In terms of the distance *ξ*=*Ut* behind the melting isotherm,
Subject to the condition *s*=*s*_{0} ahead of the melting isotherm, these equations have the solution
where *l*, the thickness scale of the melting interval, is

To gain an indication of the magnitude of *l* we take the estimates of Brown (2001) that the Melt Zone can be 15 km thick and forms in 10^{4}–10^{6} years. This gives average values of *U* between 4.76×10^{−8} m s^{−1} and 4.76×10^{−10} m s^{−1}. If we take *T*_{e}=720^{°}C (corresponding to reaction I at 1 GPa pressure), *grad* *T*=20^{°}C km^{−1} and *γ*=10^{−2} s^{−1} to 10^{−3} s^{−1} (for which there is great uncertainty) then *l* is in the range 1 mm to 1 m. Gradients in temperatures greater than 20^{°}C km^{−1} reduce these estimates marginally. If the time taken for formation of this Melt Zone is 10^{7} years as in the case of the Tuolumne intrusive suite (Coleman *et al*. 2004) then *l* is 3m. This discussion has assumed that anatexis corresponding to a particular reaction occurs at a uniquely defined *P* and *T*. However, Montel & Vielzeuf (1997) for instance show that, for reaction IV, melt is produced over a temperature range of 50^{°}C at low pressure (100 MPa) and over at least 200^{°}C at high pressure (1 GPa). As emphasized by Phillips (1991) for reactions where the melt is produced over a temperature range, the physical thickness of the melting interval will be smeared out over that range. Thus the actively melting interval is expected to be quite thick, especially at high pressures. However, for these situations, Phillips (1991) shows that smearing out the thickness of the melting interval results in the total generation of melt, and hence the melt pressure generated just above the anatectic front, being the same as when the melting interval is thin.

Note that if the isotherms rise only by conduction of heat, then the time scale for a thermal perturbation at the base of the crust to reach 15 km above the base of the crust is (15×10^{3})^{2}/10^{−6} s or 7.1 Myr; for a Melt Zone 20 km thick this time is 14 Myr. Thus, if the melt region is to reach a thickness of 15 km in 10^{4}–10^{6} years the heat must be advected upwards by moving melt. It is instructive to note that data from one of the best documented granitoid systems (the Tuolumne intrusive suite described by Coleman *et al*. (2004) and Glazner *et al*. (2004)) indicate a number of intrusions emplaced over a time period of 7 Myr. This suggests that if the Tuolumne intrusive suite was produced from a Melt Zone just over 15 km thick then the melt fluxes generated at the anatectic front were relatively small such that heat advected with melt approximately matched heat transport due to conduction. We will see that this corresponds to a Peclet number close to 1.

When the reaction involves melting with 1 m^{3} of the mineral assemblage releasing *n*_{melt} cubic metres of melt, the volumetric rate of melt generation, *Q*_{melt}, is *n*_{melt}*Q*_{s}. Thus,
And for the simple case of a well-defined *T*_{e} the rate of melt generation per unit volume is
where *s*=*s*(*ξ*) is now given by . This rate is a maximum at a distance behind the moving isotherm, and decreases to zero beyond that point. The total rate of generation of melt in this reaction zone (volume of melt per unit time per unit area of the reacting zone) is
which is independent of the reaction rate, *γ*, and the temperature gradient, grad *T*.

In what follows we assume that the flow of melt within the Melt Zone obeys Darcy’s Law. This holds so long as the melt behaves as a Newtonian viscous material. However, the experimental data (Vigneresse & Burg 2004) indicate that silicate melts behave at least as visco-elastic materials even if the viscosity is Newtonian. As the temperature decreases, non-Newtonian and ‘yield’ behaviour ultimately eventuate. The flow of non-Newtonian and visco-elastic materials through fractures and porous media is not well understood and quite recently the role of visco-plastic behaviour has been raised as a relatively new development (Ancey *et al*. 2009). A comprehensive review of older literature is given by Savins (1969). Instabilities in such materials are again poorly understood (Goddard 2003) even though it appears that the flow of these materials in fractures and porous media is characterized by instabilities such that fluid flow localizes into transient channels. In particular, behaviour similar to nonlinear Forchheimer behaviour is well known (Getachew *et al*. 1998). This means that non-Newtonian fluids flow through fractures and porous media more slowly than is expected of Newtonian fluids under identical conditions and so this needs to be taken into consideration in many of the estimates of melt fluxes made in the following since such estimates are likely to be over-estimates if the effects of melt elasticity or non-Newtonian behaviour become important (Wissler 1971).

The vertical Darcy flux of the melt, , immediately at the anatectic front is the volume of melt produced per unit time per unit horizontal area or ℜ_{melt}. Thus if *z* is the vertical distance above the Melt Zone, the vertical gradient of the total melt pressure at the top of the anatectic front is
where is the total melt pressure, *μ* is the melt viscosity, *K* is the vertical permeability, *ρ*_{melt} is the density of the melt and *g* is the acceleration due to gravity. This pressure is a constant once the isotherm has risen a distance equal to the thickness of the reaction zone and the melt pressure induced by the upward flow diffuses upwards according to a simple diffusion equation (Phillips 1991)
where *κ*_{melt} is the pressure diffusivity of the melt and *p* is the reduced pressure defined as , where *ρ*_{0} is the average density of the melt and *z* is directed vertically upwards.

Many authors have considered mechanisms that operate within the Melt Zone to drive the flow of melt. As discussed earlier, these include deformation (Brown *et al*. 1995; Brown & Rushmer 1997; Brown 2004, 2008; Rabinowicz & Vigneresse 2004; Rutter & Mecklenburg 2006), which potentially can drive melt both vertically and horizontally (table 2), and compaction (McKenzie 1984), which drives melt upwards. In addition, topographic gradients drive melt horizontally. As far as the development of plutons higher in the crust is concerned only the vertical component of melt velocity within the Melt Zone is important in contributing to the magnitude of . However, all these processes contribute to fluctuations in pressure at the anatectic front so that the above diffusion equation is replaced by
where *F* is a nonlinear function mainly defined by the deformation processes operating within the Melt Zone and λ_{c} is a parameter such as strain that describes the evolution of the deformation. This last equation is a reaction–diffusion equation and, depending on the form of *F*, may have steady or periodic solutions or much richer forms of behaviour (Cross & Hohenberg 1993). Thus steady flow or periodic pulses of magma might arise from such coupling or even chaotic behaviour. The issue deserves deeper discussion, but for the moment we consider only the simple diffusion equation. Then, at the anatectic front,
and the pressure at the anatectic front increases with time; thus
The total pressure at the anatectic front evolves as
and so is always greater than the hydrostatic melt pressure and with increase in time becomes super-lithostatic (Phillips 1991) if the initial permeability of the rocks above the anatectic front is small enough and/or *μ*ℜ_{melt} is large enough. It is easy to see, however, that if any or a combination of any of the following occur, then the pressure at the anatectic front evolves to a new value: (i) the melt-producing reaction alters thus changing the value of *n*_{melt}, (ii) the upward velocity, *U*, of the reaction isotherm changes due to the deformation processes (representing rapidly evolving perturbations) mentioned above or to additional heat input (these represent long period perturbations, and (iii) the modal concentration of the hydrous phase changes and so alters *s*_{0}. Thus the processes occurring near the anatectic front control the pressure at the reaction front and hence the pressure distribution throughout the entire granitoid system. Note also from this argument that it is the upward velocity of the melt, , that is controlled by these processes. As discussed in §6, the permeability is a dependent variable that has to evolve, by compaction or dilation, to a value that accommodates this flux given the pressure gradient that develops. This point cannot be emphasized too strongly. Traditional discussions of melt transport treat the permeability as a prescribed parameter and the melt flux as a variable that depends on the permeability. These discussions are inherited from classical ground-water hydrology where the distribution of permeability is defined by the local geology and the Darcy flux changes accordingly. In a granitoid system, the melt flux is relentlessly produced at the anatectic front and must be continuously accommodated higher in a lithostatically pressured system by the evolution of the permeability as the viscosity evolves. Concepts such as aquifers and aquitards are no longer relevant.

### (b) Permeability within the Melt Zone

We adopt the proposal made by Brown & Solar (1998*a*,*b*) that the leucosomes (defined by both lineations and foliations) have at some stage been the site of conduits for melt migration within the Melt Zone. Similar proposals are made by Collins & Sawyer (1996) that steeply plunging fold hinges act as conduits and a detailed discussion of the implications of this proposal is given by Brown (2010). We represent these conduit systems as arrays of plunging cylinders with diameter *δ* or as arrays of dipping plates with thickness *δ* and spacing *l*_{0} as shown in figure 2.

For both arrays, and for a vertical gradient in melt hydraulic potential, the permeability tensor *K*_{ij} is given by (Phillips 1991)
where *D*=*n**π**δ*^{4}/128 for an array of cylinders, with *n* being the average number of cylinders per unit area in a plane normal to the cylinders, and *D*=*δ*^{3}/12*l*_{o} for an array of plates. Thus for cylinders of diameter 10^{−2} m with *n*=10 m^{−2} and *θ*=0^{°} the vertical permeability is 2.5×10^{−9} m^{2}. If the same array plunges at 30^{°}, the vertical permeability is 6.25×10^{−10} m^{2}. For a vertical array of plates with *δ*=10^{−2} m and *l*_{0}=10^{−2} m, the vertical permeability is 8.3×10^{−6} m^{2}, whereas if the dip is 30^{°} the vertical permeability is 2.1×10^{−6} m^{2}. Thus, if one assumes that the linear and planar leucosomes act as channel ways then astoundingly large vertical permeabilities are indicated even if the lineations and foliations are gently inclined to the horizontal. Such estimates assume that all the conduits contribute to the permeability at any time; that is, there is always complete connectivity between the conduits and no dead ends exist. However, we emphasize that the permeability structure that develops throughout the entire system is not arbitrary but is controlled by and the melt pressure gradient. We will see that these simplistic models involving complete connectivity of porosity lead to permeabilities that are unrealistically high. However, one should also note that these calculations assume Darcy flow. As indicated earlier (§5*a*) nonlinear effects arising from non-Newtonian behaviour of the melt mean that could be less for these permeabilities than expected from Darcy flow.

These large estimates of permeability are comparable to those made by Petford & Koenders (1998) on the assumption that permeability in the Melt Zone is generated by brittle cracking arising from thermal stress. We return to the issue in §5*e*, where we consider advection of heat in the Melt Zone and the possibility of convective instabilities developing in the Melt Zone.

### (c) Deformation–melting coupling

Many authors have emphasized a close interdependence of deformation and melting (Brown & Solar 1998*a*,*b*). If inelastic deformation and melting are coupled so that each process influences the other then at the scale of 1–100 m the energy dissipated by the deformation may be considered as one energy source that drives the melting reaction (Hobbs *et al*. in press). Then for a power-law viscous material, the effective viscosity, *η*, is a function of the strain rate, , and can be written
where *η*_{0} is a reference viscosity equal to the viscosity at the beginning of straining and for a strain rate , *ε* is the strain and *p* is an exponent that describes strain softening (*p*<0) or hardening (*p*>0). This expression assumes that the melting reaction rate is not influenced by the strain rate. This means that the material behaves as a strain-rate softening material as well as showing strain softening or weakening. The dependence of the effective viscosity on for coupled melting behaviour at the scale of 1–100 m is identical in principle to thermal–mechanical coupling at the kilometre scale as described by Hobbs *et al*. (2008). The dependence on the inverse square of the strain rate is identical in both cases and results in shear zone localization, folds and boudins even in layered materials with low contrasts in competency where classical Biot theory predicts such structures would not form. The coupling also means that melt reactions tend to go to completion in high-strain zones before other lower strain parts of the rock mass. Figure 3 shows the results of modelling layered materials with deformation coupled to melting. The results are similar to the fabrics developed in partially melted rocks with leucosomes developing in the necks or on the margins of boudins and as zones approximately parallel to the axial planes of folds (Vernon & Paterson 2000). We propose that most of the melting fabrics observed in migmatites arise from deformation–melt coupling (see also Brown & Solar 1998*a*) so that the development of a melt ‘plumbing system’ within the Melt Zone is dependent on this coupled behaviour.

### (d) The Gurson constitutive model: ductile fracture

The growth and coalescence of small cracks or voids is well known as a mechanism of mechanical failure in rocks although most attention to this issue in geology has been in the so-called ‘brittle’ field. More general considerations of this kind of behaviour in geology are reported in Regenauer-Lieb (1999), Eichhubl & Boles (2000), Eichhubl *et al*. (2001), Eichhubl & Aydin (2003), Eichhubl (2004) and Brown (2004) and involve opening-mode ductile and brittle fractures that form by ‘void’ coalescence. The process is well known in ceramics (Ashby *et al*. 1979). In the metals literature, void growth has been studied in nonlinear (power-law) viscous materials by Budiansky *et al*. (1982), Needleman *et al*. (1995) and Subramanian *et al*. (2005). Void growth in rate-independent plastic materials (commonly labelled brittle in the geological literature) has been analysed by Tvergaard (1989) and Needleman *et al*. (1992).

Although the term ‘void’ is used in the mechanics literature the term carries the implication that the void is empty though, even in metals, the void is filled with air, or at high temperatures, a vapour in equilibrium with the metal. In order to make the terminology perhaps more palatable to geologists interested in granitoid systems we refer to these structures as ‘melt inclusions’ by analogy with the term ‘fluid inclusions’. It is understood the melt inclusions can occupy grain boundaries just as fluid inclusions can. The shape and size of these inclusions is controlled by the stress field and the competition between the surface energy of the inclusion and diffusion of melt to the inclusion (Subramanian *et al*. 2005). These inclusions presumably nucleate by mechanisms documented by Sawyer (2001) and grow by coalescence.

The constitutive behaviour is expressed by a yield function given by
where is the yield stress of the matrix material, *f* is the volume fraction of voids and *q*_{1}, *q*_{2} are constants that are normally given values of 1.25 to 2 and 1, respectively. *σ*_{e} and *σ*_{h} are the equivalent and mean stresses given, respectively, by and *σ*_{h}=*σ*_{kk}/3 where *σ*′_{ij} is the stress deviator.

Localization in these materials can occur at any angle to the stress axes and the nature, orientation and timing of localization is strongly dependent on the stress triaxiality, defined as *σ*_{h}/*σ*_{e}, and also on the detailed mechanisms by which voids grow and interact. A recent summary is given by Nahshon & Hutchinson (2008). For specimens deformed in biaxial stress states (especially if one principal stress is extensional) shear bands with void sheets form inclined to the stress axes (small or negative values of triaxiality). For large values of triaxiality (corresponding to extensile deformations) void sheets form normal to the extension axis. Ohno & Hutchinson (1984) show that localization in an extensional deformation is nucleated in regions of higher void concentration and that deformation then localizes into planar regions where the void concentration grows with continued deformation. These periods of void growth are associated with the arrest of plastic deformation in the bulk of the material which unloads to deform elastically. Void growth is greatest in the central regions of the localized band since the local effective triaxiality is greatest there. We propose, by analogy with metals and ceramics, that ultimately melt inclusions coalesce in the centralized regions of localized deformation bands to form a melt-filled fracture that may be represented as a leucosome or dyke. At high temperatures this would be expressed as a ductile dyke since it forms in the ductile regime. At low temperatures a brittle dyke forms. Thus, the mechanism of dyke formation is envisaged not by elastic fracture mechanics as proposed by many authors including Clemens & Mawer (1992) but by melt inclusion coalescence and ultimately ductile fracture (Brown 2004). The structures that form by such a process in the Melt Zone presumably resemble the dyke structures commonly observed in migmatites with intimate intermingling at their margins with the leucosomes in the adjacent migmatite (Brown 2004, fig. 3). Note that the formation of fracture systems by these processes is related to the stress triaxiality and does not necessarily conform to predictions based on either classical Andersonian or Mohr–Coulomb theory.

### (e) Advection of heat

Once a high permeability system of leucosomes is established within the melt production region, heat can be advected upward at a faster rate than simply by conduction. This is an important event in the evolution of the melt plumbing system since it enables the upward speed of the isotherms to increase substantially and hence increase and the melt pressure above the reaction front. We consider the system shown in figure 4.

Zhao *et al*. (2008) have shown that the vertical temperature distribution arising from an upward flowthrough of magnitude is given in dimensionless form by
5.1
where *P**e* is the Peclet number given by . *Pe* expresses the ratio between the transport of heat by advection relative to that by conduction. Advection exceeds conduction for *Pe*>1. *T** is the dimensionless temperature given by
5.2
and *z** is the dimensionless vertical coordinate given by *z**=*z*/*H*. *H* is the height of the system, *c*_{p} is the specific heat of the melt and λ_{e} is the effective thermal conductivity of the combined rock and melt given by λ_{e}=*Φ*λ_{melt}+(1−*Φ*)λ_{rock} where λ_{melt}, λ_{rock} are the thermal conductivities of the melt and rock, respectively, and *Φ* is the porosity; *q* is the thermal flux at the base of the subsystem of interest and *T*_{top} is the temperature at the top of the subsystem of interest. If the vertical gradient in melt pressure is lithostatic and melt-flow obeys Darcy’s Law with a vertical permeability, *K*, then is given by
In figure 5*a* we plot the dimensionless temperature profiles for various values of the Peclet number between 1 and 5 assuming a lithostatic pressure gradient for the melt pressure. *Pe* is governed by the values of *H* and and so figure 5*a* can be interpreted as showing that, for a given *H*, the temperature difference across *H* increases as increases. For values of *P**e*>5 the temperature difference across *H* increases further and the temperature rapidly approaches an isothermal distribution over most of *H*. Such high values of *Pe* have been modelled by Leitch & Weinberg (2002).

To be specific, the temperature distribution is given by (Zhao *et al*. 2008)
In the calculations that follow, the values of various parameters are assumed to be *c*_{p}=1.2 KJ kg^{−1} K^{−1} (Dingwell *et al*. 1993), *ρ*_{rock}=2700 kg m^{−3}, *ρ*_{melt}=2300 kg m^{−3} (Murase & McBirney, 1973), *μ*=10^{4.8} Pa s (Clemens 1998), λ_{rock}=3.35 W m^{−1} K^{−1} (Turcotte & Schubert 1982), λ_{melt}=1.26 W m^{−1} K^{−1} (Murase & McBirney 1973) and *q*=0.06 W m^{−2}. For simplicity, we assume throughout λ_{e}=3 W m^{−1} K^{−1} corresponding to *Φ*=0.17. Since we have assumed that the melt pressure gradient is lithostatic, *K* is fixed once is fixed. For a lithostatic melt pressure gradient and *μ*=10^{4.8} Pa s, . The temperature profiles that result assuming the anatectic front is at 700^{°}C and is 5, 10 and 15 km above the base of the crust are given in figure 5*b*. This shows that Peclet numbers in the range 1–2 are all that is required to be compatible with temperatures of 1200–1300^{°}C at the base of a 15 km thick Melt Zone with a melting temperature at its top of 700^{°}C. *Pe*=1.4 gives a temperature of 1300^{°}C at the base of a 15 km thick Melt Zone if the temperature at the anatectic front is 700^{°}C. This is probably an upper limit for what occurs in Nature. If we assume *Pe*=1 for this system then m s^{−1}; *Pe*=1.4 gives m s^{−1}. This would seem to be an upper limit for in a Melt Zone 15 km thick with a temperature at the anatectic front of 700^{°}C which for reaction I implies a pressure of 800 MPa or a depth of 30 km (Peto 1976). If we assume the temperature at the anatectic front is 600^{°}C corresponding to the wet granitic solidus at 500 MPa pressure then temperatures of 12–1300^{°}C quoted above decrease to 1100–1200^{°}C. For a lithostatic melt pressure gradient, a value of m s^{−1} gives an upper limit for *K* of 1.58×10^{−9} m^{2}. This value is significantly smaller than the estimates made in §5*b* so that if leucosomes are conduits for melt transport then the interconnection between leucosomes at any time must be small as is also emphasized by authors such as Brown (2007) and Sawyer (2001). We note in advance of §6 that these values of *Pe* are very much lower than values adopted by Leitch & Weinberg (2002), who used values of *Pe* (as defined here) in the range 50–250.

This now allows us to estimate the pressure diffusivity *κ*_{melt} within the Melt Zone. *κ*_{melt} is given by (Phillips 1991)
where *V*_{p} is the speed of sound in the melt and *μ* is the viscosity of the melt. *V*_{p} is given by where *B* and *G* are the bulk and shear moduli of the melt. If we take *B*=37 GPa and *G*=27 GPa as representative values of these elastic moduli at 650^{°}C and *ρ*_{melt}=2300 kg m^{−3} from Bagdassarov *et al*. (1992), we obtain *V*_{p}=5.6×10^{3} m s^{−1}. Taking *μ*=10^{4.8} Pa s, *K*=1.58×10^{−9} m^{2} and the geometries discussed above with *Φ*=*δ*/*l*_{o}=10^{−3}/10^{−1}=10^{−2} for planar leucosomes, we arrive at *κ*_{melt}=3.9×10^{−6} m^{2} s^{−1}, which is approximately four times larger than *κ*_{thermal} so that pressure perturbations in the Melt Zone diffuse approximately four times faster than temperature perturbations. Larger porosities result in proportionally smaller values of *κ*_{melt}.

We also need to consider whether this flow system is stable or whether conditions for convection of melt within the Melt Zone are attained. Convection in systems with upward flowthrough driven by a fluid pressure gradient greater than hydrostatic is considered by Zhao *et al*. (2008). In contrast to a system where the temperatures and pressures are fixed at both the top and base where throughflow stabilizes the flow (Jones & Persichetti 1986), a system where the basal boundary conditions comprise fixed heat and mass fluxes and the melt temperature and pressure are fixed at an upper (permeable) boundary is destabilized by throughflow so that increases in *Pe* increase the possibility of convection (fig. 6*a*, from Zhao *et al*. 2008). The Rayleigh number for such a system in which the temperature and melt pressure is fixed at the top and a heat flux, *q*, and melt pressure gradient are fixed at the base, is given by
where *β* is the volumetric thermal expansion coefficient of the melt, *ρ*_{melt}=2300 kg m^{−3}, *c*_{p} = 1.2 KJ kg^{−1} K^{−1}, *g* = 10 m s^{−2}, *q* = 0.06 W m^{−2}, *K* = 1.58 × 10^{−9} m^{2}, *H*=15 km, *μ*=10^{4.8} Pa s, and λ_{e}=3 W m^{−1} K^{−1}. Webb *et al*. (1992) report values of *β* close to 10^{−4} K^{−1} for dry sodium silicate melts. It appears that the addition of several percent H_{2}O to the melt increases this value by an order of magnitude (Ochs & Lange 1997) to values close to 10^{−3} K^{−1}. We assume for wet melts that *β*=10^{−3} K^{−1}. This gives *Ra*=2.39 and *Pe*=1.4. For values of *K* different from 1.58×10^{−9} m^{2} both *Ra* and *Pe* change proportionally.

Figure 6*a* shows that the conditions defined by *Ra*=2.39 and *Pe*=1.4 are well within the stability field for upward flowthrough systems whereas a fivefold increase in *K* to 8×10^{−9} m^{2} gives *Ra*=11.95 and *Pe*=7, which represent conditions well over the boundary of that defining convective instability and are unlikely to be reached because the system would have started to convect well before evolving to that stage. Thus one sees that the system is finely tuned; an upward melt velocity a little greater than 10^{−10} m s^{−1} for a lithostatic melt pressure gradient results in *K*>1.58×10^{−9} m^{2} for a Melt Zone 15 km high with 700^{°}C at the anatectic front and induces convection in the Melt Zone with temperatures at the base that are unrealistic. m s^{−1} under these conditions is stable and produces temperatures at the base that are at the limit of what is to be expected. It should be noted that figure 6*a* is derived analytically for conditions defined by *Pe*<6. There are no analytical solutions available for *P**e*≥6 but, as we have seen, *P**e*≥6 seems unlikely for systems of crustal dimensions.

Convective flow in upward flowthrough systems is slightly different from that developed in systems at an ambient hydrostatic pore pressure gradient in that the convective flow tends to either reinforce or compete against the upward flow generated by the non-hydrostatic pressure gradient as shown in figure 7*b* and full convective rolls develop only at large *Pe* (figure 7*c*), as discussed by Zhao *et al*. (2008). This results in localized high-velocity throughflow regions separated by lower velocity throughflow regions with elevated isotherms in the higher velocity regions. At least for systems close to the critical conditions for convective instability, the pattern is periodic in space with a horizontal wave number, *k*, given by (Zhao *et al*. 1999, 2000, 2008)
where λ is the horizontal wavelength for the upwelling regions. For *Pe*=1.4 and *H*=15 km the resulting wavelength for the upwelling regions is 19.5 km. Arguments such as these may be relevant to understanding the spacing of individual plutons within a batholith in the upper crust. Examples of flow and temperature distributions in convective systems with throughflow for different values of *Ra* and *Pe* are given in Zhao *et al*. (1999, 2000, 2008). For relatively large values of *Pe*, as represented in figure 7*c*, a convection roll forms with melt drawn down to produce colder regions and upwards to produce hotter regions. In principle, this initially convects the melt back across the path it has just travelled in *P*–*T* space but the full development of such a path assumes the melt viscosity remains constant.

### (f) Relaxation of melt pressure gradient

The pressure distribution within the Melt Zone is not static but varies discontinuously in both space and time as various reactions are initiated and as rock units with different volume concentrations of hydrous minerals (fertility) are encountered by the rising isotherms. This means that both *n*_{melt} and *s*_{0} change within the Melt Zone, resulting in fluctuations in ℜ_{melt} at the top of the Melt Zone and in the pressure gradient within the Melt Zone. From time to time, especially when thick units of relatively infertile compositions (low *s*_{0}) are encountered, the supply of melt decreases so that the vertical gradient in drops below lithostatic. When this occurs, although the new mean total pressure in any particular compartment must be lithostatic, the pressure gradient relaxes to hydrostatic on the time scale of (*L*^{2}/*κ*_{melt}). This produces compressive stresses at the base of the compartment and tensile stresses at the top as discussed by Zhao *et al*. (2008). The maximum height, *H*_{critical}, of a body of rock that can support a hydrostatic pressure gradient is
where and are the tensile and compressive strengths of the melt-bearing rock mass. Thus, if MPa say, then *H*_{critical}=12.5 km and the timescale for relaxation of the melt pressure gradient to hydrostatic in this compartment is approximately 1.1 Myr assuming *κ*_{melt}=3.9×10^{−6} m^{2} s^{−1}. If the height of the Melt Zone is greater than *H*_{critical} then the base of the column collapses and the top localizes in melt inclusion sheets resulting in compaction of the rock mass at the base and expulsion of melt at the top of the compartment. Thus the Melt Zone compacts from its base upwards and releases melt at the top. It is to be emphasized that this compaction process is the result of melt migration and not the cause, as is sometimes proposed in melt systems (McKenzie 1984). These processes of course are especially pertinent as the system is closing down.

## 6. The mechanics of the Transport and Emplacement Zones

A model for the Transport and Emplacement Zones has been presented by Leitch & Weinberg (2002). We summarize that paper here to indicate the differences between that model and the one presented here and to highlight the importance of the parameter *Pe*. Leitch and Weinberg model the Transport and Emplacement Zones as their ‘free rider layer’ within which the melt is supra-solidus but loses latent heat as it rises. Their model results in most of the free rider zone being isothermal, with a rapid decrease in temperature towards the top before blending into the ambient geothermal gradient above the layer. The magnitude of the rapid decrease in temperature increases with time and with the magnitude of *Pe*. In their model, thicknesses of the free rider zone can reach approximately 40 km within 0.5 Myr.

It is important to note that the Peclet number is defined by Leitch and Weinberg in terms of the physical velocity of the melt, *V* , instead of the Darcy flux, , as in this paper. The two are connected by so that a Peclet number of 1000 quoted by Leitch and Weinberg together with a porosity of 0.05 is equivalent to *Pe*=50 in this paper. We have seen that values of *Pe* an order of magnitude smaller than 50 are already large numbers in a crustal melt system. Figure 5*a* shows that, for *Pe*=5, an isothermal lower region of a system develops. Most of the effects (including the very large thicknesses of the free loader zone) and the rapid rates for their formation modelled by Leitch and Weinberg arise from the very large values of *Pe* assumed (*Pe*=12.5–200, using in defining *Pe*). These issues are important to grasp because the essence of the principle of maximum entropy production rate that will be discussed in §7 is that systems driven from equilibrium by heat input cannot grow indefinitely and new processes arise, driven by the available energy sources, to compete with existing processes to prevent the development of such large Peclet numbers that correspond to dissipation processes that are physically unattainable. We present below a different model of the Transport/Emplacement Zones where it is convenient to distinguish two subzones that merge into each other with increasing height. In doing so, some of the important processes that control the nature of the melt plumbing system are highlighted. The essence of the model presented here is that the evolution of the Transport and Emplacement Zones is controlled by a well-constrained flux of melt from the anatectic front and this results in values for *Pe* much lower than those assumed by Leitch & Weinberg (2002).

We assume that the constitutive behaviour of the rocks constituting the Transport Zone can be described by the Gurson–Tvergaard yield function even though the rock mass itself deforms as a rate-dependent (viscous) power-law material at its base near the anatectic front and as a rate-independent (brittle) material towards its top. Then localization and its geometry are dependent on the amount of strain, the triaxiality of the far field stress and, most importantly at the base of the zone, on the local concentration of melt inclusions (Ohno & Hutchinson 1984). This means that in compressional regimes the localized melt inclusion sheets will tend to occur as shear zones inclined to the principal compressive stress whereas in extensional regimes they will tend to be vertical zones normal to the minimum principal stress. In all cases though, localization will develop preferentially in regions at the anatectic front with higher concentrations of melt arising from bodies of rock with higher *s*_{0} or from upward convective welling of melt. Having established a system of localized melt channels the pressure distribution at the anatectic front diffuses upwards through a distance, *L*, on a time scale given by (*L*^{2}/*κ*_{melt}) as described by Phillips (1991). As discussed above, for a Transport Zone 10 km high this time scale is of the order of 0.7 Myr. This sets a limit for the minimum time that is available to transport a first batch of melt from the anatectic front to the top of a 10 km thick Transport Zone. For a Transport Zone 20 km thick this time is 2.9 Myr. Note that the time scale of 0.5 Myr associated with the growth of a free rider layer 40 km thick modelled by Leitch & Weinberg (2002) implies a relatively large melt diffusivity of 10^{−4} m^{2} s^{−1} as opposed to a value of 3.9×10^{−6} m^{2} s^{−1} proposed here.

Melt can continue to move upwards only if the pressure is sufficient to hold the walls of the localized zones open and the gradient in hydraulic potential is sufficient to drive the melt upwards. We assume these two conditions are met if the melt pressure is lithostatic throughout although it is probable that it is super-lithostatic for some distance above the anatectic front (Phillips 1991). An additional requirement is that the melt temperature is greater than the solidus of the melt. This means that heat has to be advected upwards with the melt such that the temperature at the base of the combined Transport/Emplacement Zones is equal to the melting temperature whilst the temperature at the top is equal to that of the melt solidus. As long as the melt transport is in relatively closely spaced fractures and the melt is in local thermal equilibrium with the host rocks, this temperature difference, Δ*T*, can be obtained from equations (5.1) and (5.2) for *z*=*H*
5.3
where *q* is the heat flux into the system.

Δ*T* is plotted against *H* for various values of and *q*=0.06 W m^{−2} in figure 8. The maximum value of Δ*T* in the system is perhaps 150^{°}C (at a maximum) so that figure 8 predicts a maximum vertical height for the Transport/Emplacement Zones of about 7.5 km for m s^{−1}, and 9 km for m s^{−1}.

In the model presented here the temperature gradually decreases with height in the system instead of remaining isothermal with a temperature discontinuity at the top as arises for high values of *Pe* (Leitch & Weinberg 2002). This means that the melt viscosity decreases as one proceeds higher in the system and according to Baker (1998) undergoes a rapid decrease at temperatures governed by the H_{2}O content. Darcy’s Law demands that, in order to keep the melt flux at , the permeability of the system must increase by slightly more than the same amount that the viscosity decreases. We have seen that for a dyke-like feature (two parallel plates) the permeability is given by *K*=*δ*^{3}/12*l*_{o} and, hence, (∂*K*/∂*δ*)_{lo}=*δ*^{2}/4*l*_{o}. Three processes now operate if there is a population of conduits with different *δ*, as follows: (i) The cubic dependence of permeability on thickness of the conduit means that melt is focused preferentially into the thickest conduits (Zhao *et al*. 2008). The high-permeability conduits attract most of the melt and the low-permeability conduits are starved of melt and collapse in thickness. The result of such collapse processes would be thin granite ‘veins’ that are the remnants of collapsed conduits (Bons *et al*. 2008). This process of melt capture by high-permeability conduits and starving of low-permeability conduits is similar in principle to that illustrated by Brown (2010, fig. 9) except that in a lithostatically pressured system the ‘starved’ conduits must collapse. (ii) When an increase in permeability is required to accommodate at an increased melt viscosity, the dependence of on the square of the thickness means that an increase in permeability can be achieved readily for conduits where the initial thickness is small compared with the spacing. However if the initial thickness is large with respect to the spacing, or comparable to the spacing, then an increase in thickness means that many adjacent conduits will amalgamate. As an example, conduits 10 cm wide and spaced at 1m result in a permeability of 8.3×10^{−5} m^{2}. This is also achieved by conduits 3.16 cm thick spaced at 3.16 cm. A 10-fold increase in permeability can be achieved by increasing the thickness of the 10 cm wide conduits to 20 cm with a spacing of 80 cm. This same new permeability can be achieved by the 3.16 cm wide fracture increasing in thickness to 4 cm spaced at 6 mm. With any further increase in permeability the 4 cm wide conduits must amalgamate with adjacent conduits in an attempt to achieve the desired result. Hence conduits initially closely spaced grow fast and some will coalesce, and initially widely spaced conduits remain widely spaced. The result of this coalescence process would be thin screens of ‘host rock’ between, and embedded in, what are now the thicker granite-filled conduits. Examples of such fabrics might be figs. 2–4 in Brown & Solar (1998*b*). (iii) The strain associated with widening of a single conduit to accommodate an imposed melt flux at a given viscosity is independent of its current thickness, but because more thin conduits are required than thick conduits to produce a given permeability, the work done in widening a thick conduit array is less than that required for an array of thin conduits.

This means there is a competition between these three processes: the work done to increase the thickness in an array of thick conduits is less than that for an array of thin ones for a fixed permeability. However, the thin ones are more likely to coalesce than thick ones for a given permeability. Meanwhile the individual high-permeability conduits at any particular time attract most of the melt flow whilst the others are starved of melt and collapse. The evolution of this system is dynamic and needs further analysis to make definite predictions. However, it seems that the ultimate result is a coarsening of the overall structure, mainly by coalescence of conduits, which means that an array of thick, widely spaced conduits develops. The melt transport mechanism switches to thicker conduits that enable continuity of melt flux to be maintained. The ascent process then is that described by Clemens & Mawer (1992) and Petford *et al*. (1993, 1994) in that the melt is no longer in local thermal equilibrium with the host rocks and loses heat to its surroundings as it ascends and progressively crystallizes. The control on a switch in transport mechanisms is a decrease in melt viscosity. The melt can continue its upward transport only as long as it can advect heat and retain its latent heat of fusion.

However, even another process may influence the plumbing system as the magma cools. If the magma viscosity becomes non-Newtonian and particularly if some form of ‘yield’ behaviour ensues then the thickness of an established conduit must expand even further to accommodate the melt flux. This represents a discontinuity in the growth of the thickness of the conduits as the conditions for a change in melt constitutive behaviour are encountered. In summary, the decrease in viscosity and changes in constitutive behaviour of the magma as it rises and cools result firstly in a selection process operating as the viscosity first starts to increase and then a discontinuous increase in thickness of surviving conduits as the mechanical properties of the melt change character. This is illustrated in figure 9.

Thus the Melt Transport Zone may consist of two subzones that blend into each other: a lower one where the transport network is closely spaced and local thermal equilibrium with the host rocks is attained so that, on a regional scale, isotherms continue to be elevated. Here the temperature is supra-solidus for the melt but sub-solidus for the host rocks. Above that subzone the transport fracture network is widely spaced and the melt loses heat to its surroundings with elevated isotherms only in the vicinity of fractures. The variety of plutons described by Brown & Solar (1998*b*) from a single terrain seem to have roots in both of these subzones.

Horizontal flow induced by topographic gradients can use mechanical weaknesses in horizontal bedding parallel foliation in the lower subzone in the Transport Zone. Given the horizontal melt pressure gradient quoted in table 2, a melt viscosity of 10^{4.8} Pa s and a melt flux of 10^{−10} m s^{−1}, we arrive at a permeability of 4.6×10^{−9} m^{2}. It is interesting to note that this permeability can be achieved by melt channels 8 mm thick separated by 10 mm of host rock. Thus *lit par lit* fabrics may represent the development of such permeability networks. The Peclet number for horizontal flow over 20 km (the dimension quoted for near horizontal flow in the Himalayas by Leitch & Weinberg (2002)) is between one and two times that for vertical flow driven by buoyancy. This means that the horizontal flow is capable of advecting considerable heat horizontally and in regions of substantial topographic gradients the ‘thermal plume’ arising from the Melt Zone is inclined upwards in the direction of horizontal flow in a manner implied by fig. 2 of Brown (2010).

The theory behind the discussion presented above is now presented. Once an upward melt flow regime (comprising an upward flow and associated melt pressure gradient and permeability) is established just above the anatectic front the permeability is self-adjusting throughout the Transport/Emplacement Zones and depends on the upward melt flux (measured in cubic metres per square metre per second). As discussed by Zhao *et al*. (2008, pp. 10–13), the permeability established just above the anatectic front acts as a pressure valve for the system above it. If this permeability is *K*_{1} and the melt encounters a region above where the permeability is *K*_{2} then, since for mass continuity the upward velocity must be the same in both layers,
where we assume that the melt has the same density in both layers; *μ*_{1} and *μ*_{2} are the respective melt viscosities. It follows that if *μ*_{1}=*μ*_{2}
so that approaches the melt hydrostatic pressure gradient, *ρ*_{melt}*g*, as *K*_{2} becomes larger than *K*_{1}. For instance if *K*_{2}=10*K*_{1} and is lithostatic then is only 1.017 times the hydrostatic melt pressure gradient so the effect is dramatic.

Following the discussion in §5*f* the more permeable unit must compact from its base upwards until the permeability reaches *K*_{1}. Conversely, if *K*_{2}<*K*_{1}, the melt pressure must in principle be greater than lithostatic but now this pressure is capable of deforming the walls of the conduit so that the permeability increases to *K*_{1}. Thus the permeability within the transport zone is self-adjusting and tends towards a permeability distribution that can accommodate the flux established at the anatectic front.

On the other hand, if the melt is rising in a zone that has a constant permeability, *K*_{1}, and the viscosity changes from *μ*_{1} to *μ*_{2} higher in the system because the temperature and/or the volatile/crystal content of the melt changes, then
so that a 10-fold increase in viscosity arising from a decrease in temperature for a lithostatic value of () results in a melt pressure gradient that is more than 10 times the lithostatic gradient, which means that the melt pressure at the top of this part of the transport zone increases dramatically. This is presumably capable of initiating classical brittle fractures in the region above in a manner discussed by Clemens & Mawer (1992) if the ambient temperature of the crust is low enough but can also be accommodated by an increase in permeability by expansion of conduits as discussed earlier. This increase in viscosity due to cooling of the melt is the fundamental step in initiating the transition from flow in a closely spaced network initiated by ductile fracture to flow in widely spaced dykes as illustrated in figure 9.

At any stage during the upward transport of melt, it is possible for the melt to move horizontally into a low-strength zone to form a sill-like structure. Within this sill, although the melt pressure is lithostatic, the melt pressure gradient quickly relaxes to hydrostatic. So long as the melt constitutive behaviour is that of an elastic–viscous material the sill must collapse from its base upwards until most of the melt is expelled. The time scale for this depends on the viscosity, but, for viscosities of 10^{5} Pa s, this time is measured in days at most since for stress differences of only 1 Pa the strain rates are very fast at about 10^{−5} s^{−1} assuming Newtonian viscosity.

In order for such a sill-like structure to remain a finite thickness until it solidifies, the mechanical properties of the magma must involve yield behaviour so that the magma can support finite vertical stresses over extended periods of time. Several authors have proposed that the melt–crystal mixture behaves as a Bingham material (Petford 2003) but this alone is not sufficient to allow finite thickness horizontal intrusions to develop.

There is an important distinction between a ‘liquid’ and a ‘solid’ to be made although recent developments in non-Newtonian fluids are making this distinction more difficult to define (Ancey *et al*. 2009). The term ‘yield’ is used in two different ways depending on whether the material is a liquid or a solid. In the literature on elastic–viscous liquids, the term yield is used to describe a rapid increase in viscosity with increase in shear stress or strain rate (Barnes 1999). This behaviour is commonly modelled as a Bingham material although other nonlinear behaviour is possible. This means that the liquid continues to deform in a permanent manner when stressed both above and below the yield point, where the viscosity is large, although the strain rates are dramatically different. Silicate melts show this kind of behaviour as discussed by Petford (2003) and Vigneresse & Burg (2004). We refer to the stress where this kind of liquid-type yield becomes apparent as .

Solids show a transition between stress states where the response is elastic (or, strictly, hyper-elastic as defined by Houlsby & Puzrin (2006) so that the strain is related to the stress through a potential) to stress states where the response is a permanent strain. The transition is referred to as the yield stress and, once that stress is attained, the material can support a stress at least equal to the yield stress indefinitely unless some elastic weakening mechanism operates such as damage accumulation. As a silicate melt begins to crystallize and become a magma (melt + crystals), a stage is eventually reached where the material can support a noticeably higher stress. There is much discussion in the literature as to whether this corresponds to a critical proportion of residual melt or some other cause. Presumably as a melt solidifies to become a solid there is a transition from an elastic–viscous material (a liquid) to an elastic–plastic–viscous material (a solid) with melt inclusions (a Gurson–Tvergaard material) and then to an elastic–plastic–viscous material with no melt inclusions. As far as the development of plutons is concerned, the important transition is that from an elastic–viscous material with no plastic yield to an elastic–plastic material with a finite yield strength that we will call . Whether this transition is equivalent to the development of a critical melt fraction is not known. For a pluton to develop solely from the development of a liquid yield, , would require that the viscosity increases by at least eight orders of magnitude from 10^{5} Pa s to greater than 10^{13} Pa s so that the magma-filled space occupied by the pluton can remain open before it solidifies. A simple increase in viscosity of this magnitude seems unlikely because it would severely decrease the flow rates in supply conduits during emplacement unless very large effective permeabilities can develop.

The argument presented in §5*f* involving *H*_{crit} is applicable here. If the mechanical behaviour of the crystallizing magma is characterized by the values of the magma tensile strength, *σ*_{tensile} (which is presumably small), and yield stress, , then the critical vertical thickness of a horizontal intrusion is given by
Hence a pluton 1 km thick implies, at a minimum, a value of 4 MPa for with proportionally larger values for thicker single intrusive event plutons. A minimum value of 20 MPa for is needed to preserve a pluton 5 km thick in a single intrusive event; this relatively small value for does not seem unreasonable (Rutter & Neumann 1995; Brown & Rushmer 1997).

Thus the Emplacement Zone represents a region at the top of the system where the magma temperature is decreasing towards the melt solidus and the magma is crystallizing. The constitutive properties of the melt–crystal mixture change from elastic–viscous to elastic–plastic–viscous as the temperature decreases. The total system is still characterized by an upward melt flux, , controlled by the processes at the anatectic front, and the melt pressure remains lithostatic (or above) in the vertical vents. However, any horizontal excursion of the melt facilitated by a low-strength horizontal zone produces a region where the mean melt pressure is lithostatic but the melt pressure gradient rapidly decays to hydrostatic. This region can remain open up to a critical thickness governed by the value of for the magma. This region can expand to this critical thickness driven by the flux, .

Even though the melt pressure gradient within the pluton is hydrostatic, the mean pressure within this pluton is lithostatic, and so the shape of the pluton is controlled by the regional stress field with the pluton growing fastest in the direction of minimum principal stress. This is in agreement with observations of pluton shapes in three dimensions which appear to be controlled by the kinematics of the regional crustal deformation (Vigneresse 1995*a*,*b*; Brown & Solar 1998*a*,*b*).

## 7. The thermodynamics of the granitoid system

We have seen that the behaviour of the complete granitoid system is controlled by the conditions at the anatectic front. This is where is established, controlled by values of *s*_{0}, *U* and *n*_{melt}. The height of the transport/emplacement region is then controlled by with a relation between and *H*, expressed as the Peclet number, that defines the conditions under which the isotherm corresponding to the melt solidus is elevated to higher levels in the crust. Thus in considering the controls on entropy production rates in the system it is necessary to consider only the Melt Zone with a particular emphasis on the anatectic front. The base of the melt production part has a temperature of *T*_{BASE} at any particular time and a temperature at its top given by *T*_{REACTION}, where the term reaction refers to the melting reaction at the anatectic front. The base of the system receives a heat flux from the mantle, *q*_{IN}, and a heat flux, *q*_{OUT}, leaves this part of the system at its top. There is also dissipation of energy throughout the crust by inelastic deformation, melting and fluid flow; latent heat is absorbed in the melt region. The heat input at the base of the system, plus heat generated by internal radioactive decay, plus the deformation in the melt region drives the melting process and supplies the energy to move the melt higher into the crust. This part of the granitoid system is shown in figure 10.

If the Earth is static with no deformation or melting, then the heat introduced at the base is conducted to the surface and the entropy production rate, , is a minimum (Kondepudi & Prigogine 1998).

Following an argument presented by Ozawa *et al*. (2003), the entropy production rate, , for the system shown in figure 10 is
where *W* is the work done inside the Melt Zone by processes such as deformation and melt production and flow and *Φ* is the area of the base of the Melt Zone. Thus also increases with *Pe* in an identical manner to Δ*T* but can be reduced or reinforced by processes that operate within the system through the internal work term, (*W*/*T*_{REACTION}). We have not included the entropy production rate arising from heat generated by radioactive decay in this equation since it simply adds a constant to . Δ*T* increases with *Pe* as given by equation (6.1) and as shown in figure 8. For a constant *H*, *Pe* is controlled by , which is controlled by conditions at the anatectic front.

One could conclude that increases indefinitely as Δ*T* increases but a limit is set by the competing internal work terms that arise from instabilities that grow from increases in heat supplied to the system. In particular, this involves the onset of convective instabilities in the Melt Zone (figures 5*b*, 11*a*,*b*). Two closely related processes operate to decrease the rate of entropy production.

*Convective instabilities.*The initial onset of convective instabilities (figure 7*b*) reduces Δ*T*locally (Zhao*et al*. 1999, 2000, 2008), resulting in decreasing locally. However, a far more potent mechanism for decreasing comes into play as*Pe*increases further and especially if melt begins to be drawn down from the anatectic front as shown in figure 7*c*.*Retrograde mineral reactions and viscosity increases.*Figure 7*b*,*c*shows that, for Peclet numbers above the critical Rayleigh number, melt at the anatectic front starts to be taken on a path that takes the melt back across*T*_{REACTION}for the reaction defining the anatectic front. Such a path cannot continue very far because it results immediately in retrograde reactions.

If a retrograde reaction involving the production of hydrous phases occurs, then the concentration of H_{2}O in the melt decreases at constant temperature (Baker 1998; Dingwell 1999). This results in an increase in viscosity of the melt of several orders of magnitude as well as producing significant changes in *ρ*_{melt} and *β*; for instance, the analysis of Dingwell (1999) shows that at 730^{°}C the decrease in H_{2}O content of the melt from 4 per cent by weight to 0.5 per cent by weight results in a viscosity increase from 10^{5} Pa s to 10^{8} Pa s. The system suddenly comes to a halt and drops towards minimum values associated with heat conduction processes until a new permeability structure can be established. Retrograde reactions have been described, particularly at the anatectic front (Kriegsman 2001; Brown 2002), consisting, for instance, of an anhydrous assemblage plus *m**e**l**t*→*h**y**d**r**a**t**e*-bearing assemblage, as recorded, for example, by the common occurrence of retrograde muscovite overprinting sillimanite. This may be an expression of effects arising from convection of melt in systems that have just passed the peak of entropy production rate. Once a new permeability structure is established, the isotherms rise again with the opportunity for the newly deposited hydrous phase to undergo dehydration, thus producing melt yet again and ultimately re-initiating the convective instability. Thus a see-saw process develops at the anatectic front consisting of repetitions of
as shown in figure 11*b*. This process sets the upper limit to the thickness of the Melt Zone and repeatedly returns the system to a state of maximum entropy production rate with new pulses of melt liberated from the anatectic front. This see-saw effect continues as long as heat is supplied to the system and fertile bodies of rock exist to supply melt within the system.

The Peclet number involves the product, , and so this introduces a difficulty in defining an optimal value of *H* for a given and vice versa. In principle, there seems no reason why *H* cannot increase indefinitely as increases and it is only through sensitivity analyses such as presented in figure 5*b* that one discovers physical limitations to the value of . An extra constraint is needed to define an optimal height of the Melt Zone and this is given by assuming that the system maximizes the rate of entropy production subject to the constraints imposed on it. This is also the role that this principle plays in other systems such as the Earth’s climate system (Paltridge 2001), fluid flow in deforming porous media (Coussy 2004) and inelastic deformation of rate-insensitive materials (Houlsby & Puzrin 2006).

If we identify the optimal conditions for the granitoid system as that state that maximizes , as shown in figures 6*b* and 11, then this state is defined by *Ra*=*Ra*_{critical} from figure 6*a*. The curve shown in figure 6*a* has a very complicated analytical expression (Zhao *et al*. 2008) and so we have fitted a simple curve given by

This means that
Hence, given a particular , which is established at the anatectic front by the physical nature of the particular granitoid system (*s*_{0}, *n*_{melt} and *U*), the optimal height that maximizes and hence optimizes the system can be established. This is done in figure 12 for m s^{−1}, which gives *H*_{optimal}=21 km.

## 8. Concluding remarks

A model for granitoid systems has been presented that shows that such systems are optimized when a state of maximum entropy production rate is achieved. The system and its evolution are completely controlled throughout its life by processes that operate at the anatectic front. This includes, in particular, the generation of a melt flux, , that controls the permeability distribution in the system, the height of the transport/emplacement region and the size of the plutons. The important steps in this process are as follows:

The initiation of melting, which we propose is strongly coupled to the deformation and results in the structures (shear zones, folds, boudins, leucosomes) that are characteristic of anatectic terrains. The mechanical behaviour is described by a Gurson–Tvergaard constitutive law where melt-filled inclusions form and coalesce into regions of localized deformation to produce leucosomes and dykes.

The initiation of advection of heat upwards with the moving melt in deformation-induced channelways. This is a fundamental step that allows large volumes of melt to be produced rapidly and transported. It marks the departure from a state where the rate of entropy production is a minimum to a route where the entropy production rate can be maximized.

The migration of melt, driven in any direction by melt hydraulic potential gradients that arise from deformation. These gradients can be very significant over short distances (1–10 km) and can compete with buoyancy and topographically induced gradients on the scale of the granitoid system. Horizontal flow equal in magnitude to vertical flow induced by buoyancy can also occur driven by topographic gradients at the surface of the Earth.

The initiation of a size selection process in the melt transport conduits. This process results in the conduit system progressively coarsening upwards leaving collapsed conduits behind and progressively widening dykes high in the system.

The attainment of critical values of the Rayleigh and Peclet numbers that define the initiation of convective instabilities in the Melt Zone. This represents a state of maximum entropy production rate for the system as a whole and governs the maximum possible size of the system and the maximum attainable melt flux. From this step onwards processes that are driven by increases in energy input to the system, such as retrograde mineral reactions that increase the melt viscosity, operate to decrease the entropy production rate.

The transition from elastic–viscous behaviour for the melt to elastic–plastic–viscous behaviour for the magma as the temperature corresponding to the melt solidus is approached at the top of the system. This marks the fundamental step that allows finite thickness plutons to develop and defines the level in the crust where emplacement occurs.

The attainment of a Peclet number (i.e. an energy level) where the melt at the anatectic front is drawn down so that it undergoes retrograde reactions to produce hydrous phases. This dramatically decreases the viscosity of the melt and leads to lower entropy production rates, although the system may see-saw through this stage, producing pulses of melt, depending on the availability of heat and melt supply for the system.

Our analysis indicates that a melt flux at the anatectic front of the order of 10^{−10} m s^{−1} is close to the maximum possible given the physical parameters of the system. This corresponds to a maximum thickness for a Melt Zone of 21 km and is sufficient to produce plutons just over 3 km thick in a single intrusive event over a period of 10^{6} years. Smaller values of the melt flux at the anatectic front result in longer time scales for the existence of the granitoid system but a limit is set by melt fluxes which correspond to *Pe*=1. This situation corresponds to heating of the crust at the same rate as thermal conduction.

Granitoid systems are to be viewed as dynamic systems that continuously adjust to changing conditions at the anatectic front. The permeability throughout the complete system dynamically adjusts to these conditions. The system is finely tuned and its optimal configuration is that which maximizes the entropy production rate for the conditions generated at the anatectic front.

## Acknowledgements

We thank Mike Brown, Ed Sawyer, Ron Vernon and Mike Williams for comments that greatly improved the paper.

## Footnotes

One contribution of 17 to a Theme Issue ‘Patterns in our planet: applications of multi-scale non-equilibrium thermodynamics to Earth-system science’.

- © 2010 The Royal Society