Stepwise magma migration and accumulation processes and their effect on extracted melt chemistry

Numerical and analogue models suggest that melt production, its segregation from the solid matrix and subsequent transport and accumulation are highly dynamic and stepwise processes exhibiting scale invariant patterns in both time and length scales, which is characteristic of self-organized critical systems. This phenomenon is also observed in migmatites at several localities, where the leucosome thickness statistics obey power laws. Stepwise melt transport and deformation-enhanced melt mobility affect melt production dynamics by determining the distribution of extracted melt batch sizes and residence times of melt pockets within the host rock, which in turn would influence the geochemistry of extracted melts. We introduce a numerical approach, which enables qualitative and quantitative assessment of the effects of stress-induced melt migration and accumulation on the chemistry of partial melts. The model suggests that apart from different sources and melting percentages, deformation can be an important factor in producing geochemical variations within and between intrusive/extrusive complexes.


INTRODUCTION
According to widely accepted theory, partial melting of source rocks is the main mechanism of magma formation in the Earth's crust.The magma generation processes cover several orders of magnitude on the length scale, from the micrometre-size melt droplets at the onset of melting to the emplacement of kilometre-size granite bodies in the upper crust.
Magma generation starts with the extraction of melt from the solid matrix, a process called melt segregation (e.g. Brown 1994;Sawyer 2001) where grain-scale microscopical melt pockets residing at grain boundary junctions (e.g.Bulau et al. 1979;Jurewicz & Watson 1984;Laporte & Watson 1995) are concentrated into leucosomes.This is done mainly by coaxial deformation (i.e.pure shear), which has the effect of forcing the melt into melt-filled low-stress regions oriented perpendicular or at a high angle to the direction of compression (Bagdassarov et al. 1996;Vigneresse et al. 1996;Vigneresse & Burg 2000).The mobility of melt and its extraction from the rock are enhanced by non-coaxial tectonic forces resulting in gradients in the stress field (Vigneresse & Burg 2000;Marchildon & Brown 2001;Bons et al. 2004), which allow the melt to migrate over larger distances, accumulate and finally ascend into upper levels of the crust.Tectonic forces and deformation thus play a crucial role in driving melt segregation and transport processes.
The exact mechanism of how magma is extracted from the source and transported to upper levels of the crust without cooling and interacting with ambient crustal rocks is still a matter of debate (e.g.Petford & Koenders 1998;Brown 2004;Bons et al. in press).Many current theories of magma transport and accumulation assume a steady-state conductive magma flow network at the outcrop scale (e.g.Kruhl 1994;Brown et al. 1999;Tanner 1999;Weinberg 1999;Marchildon & Brown 2003), where a branching interconnected network of stromatic leucosomes feeds larger dykes which serve as vertical ascent channels.Mesoscale permeability of the magmatic system is achieved when the partial melting percentage exceeds the melt escape threshold -15-20 vol% of melt (Vigneresse et al. 1996), breaks the cohesion between mineral grains, and allows the melt to escape from the source.Therefore, relatively high bulk melt fractions are needed to initiate melt migration.
An alternative mechanism has been proposed (e.g.Bons & van Milligen 2001;Bons et al. 2001;Bons et al. 2004) by which magma transport and extraction take place at low bulk melt fractions and no interconnected melt network is needed.According to this concept, melt segregation and transport occur discontinuously in melt batches; magma transport over large distances is accomplished in the form of hydrofractures (Bons 2001) or self-propagating dykes (Clemens & Mawer 1992).It has been shown that a fluid-filled fracture can move by opening at one end of the crack and simultaneous closure at the other end (Weertman 1971;Secor & Pollard 1975;Takada 1990).Hydrofracture mobility can occur when the effective normal stress gradient along the fracture exceeds a threshold value.The effective normal stress gradient can be due to the difference in density between fluid/melt and wall-rock, and/or due to a stress gradient in the wall-rock.Using experimentally established and theoretically calculated values for effective normal stress and fracture toughness, it can be calculated that a water-filled fracture becomes unstable when it is longer than 5-100 m, while a vertical melt-filled fracture can reach lengths of 15-300 m (Secor & Pollard 1975;Bons & van Milligen 2001).It can also be shown that if the stress gradient is in the order of 1-10 MPa/m, fractures as short as 10 cm may become unstable and therefore mobile (Bons & van Milligen 2001).Hence, decimetrescale melt-filled fractures can be mobile at tectonically reasonable stress gradients.Natural examples of such features can be easily found in migmatite terranes where small veins are located at their formation site, while larger veins have moved away and merged with others.When melt veins merge, a larger new melt vein is formed that needs an even smaller stress gradient to move further.An important consequence is that the new melt vein can move more rapidly and more frequently per unit of time, leading to a greater likelihood of merging with other melt veins (Bons et al. 2004).The movement of the melt is thus envisaged as a stepwise movement of different melt batches, rather than continuous flow of melt through connected conduits.
Analogue experiments (Bons & van Milligen 2001;Urtson & Soesoo 2007) with sand and carbon dioxide as analogues of crustal rock and melt, respectively, support the model of discontinuous and stepwise extraction of a liquid phase from the solid matrix.These experiments as well as numerical models (e.g.Vigneresse & Burg 2000;Bons et al. 2004) suggest that a stepwise melt transportation system may quickly evolve into a selforganized critical state (Bak et al. 1988;Bak 1996) which exhibits fractal patterns and scale-invariant behaviour in both time and space-1/f (Bak et al. 1987) volume fluctuations and power law (e.g.Bonnet et al. 2001) batch size statistics.The dynamics of such a self-organized critical system cannot be described by linear transport equations but follow more complex rules (Bak 1996;Bons & van Milligen 2001).
It has been shown (Bak 1996;Bons & van Milligen 2001) that a self-organized critical transport system, once it has achieved the critical state, is capable of transferring any additional mass without modifying itself and that the structure of the system is maintained regardless of significant volume extraction.The size of leucosomes and their spatial distribution in migmatites can thus still carry information about the dynamic state of the system while it was still partially molten, although probably large volumes of magma have passed through the rock, the traces of extracted melt pools are not retained, and the volume of former melt present in the system is difficult to detect (e.g.Sawyer 2001;Bons et al. 2008).
It is easy to see on a conceptual level that the stepwise melt migration and the merging history of melt batches will influence the chemistry of the resulting melts.In order to understand this process throughout the magmatic evolution and quantify the melt batch chemical history, a numerical model presented in this paper can be of help.

ANALOGUE MODELLING OF MAGMA TRANSPORT AND ACCUMULATION
Analogue experiments have been performed to model liquid phase (or melt) segregation from the solid matrix (host rock) and its subsequent transport (Bons & van Milligen 2001;Urtson & Soesoo 2007).The setup of the experiment in a flat glass tank (Fig. 1a) and chosen analogue materials, wet sand and carbon dioxide produced by fermentation of the yeast, as analogues of crustal rocks and partial melt, respectively, allow realtime observation and recording of the system evolution and later statistical analysis of the system state at any selected timestep.It should be stressed that such experimental systems are analogues to real partially molten systems only at conceptual level -due to higher contrast between density and viscosity of liquid and solid phases, greater relative volume of generated 'melt', etc.Nevertheless, these experiments provide valuable information about the fundamental aspects of the dynamics of melt generation and transport and allow parallels to be drawn with real magmatic systems.Here we give a short overview of the modelling performed previously and discuss briefly the key aspects of these experiments.Detailed description of the experimental setup and data analysis methods can be found in Bons & van Milligen (2001) and Urtson & Soesoo (2007).
The experimental system exhibits different transport modes in order to adapt to increasing gas flux during the experiment.At a low gas production rate the gas percolates through the interconnected pore space between sand grains.The liquid percolation threshold (e.g.Vigneresse et al. 1996) is exceeded at an early stage of the experiment, determined by the packing of the sand with about 30% open porosity between the grains and the liquid-gas wetting angle.At this stage the continuous production of gas is inferred only from the steady rise of the liquid level in the tank and the escape of small gas bubbles on the surface; no apparent modifications appear in the sand column.The system is thus capable of accommodating all the produced gas in the pore space and transferring the low flux of excess gas.To transfer more gas at increasing rates of gas production, the system switches to the ballistical transport mode (term coined by Bons & van Milligen 2001) and open cracks and gas batches appear.In this mode the transport of the gas is intermittent and stepwise and occurs as discrete events by merging of gas batches, draining the gas into adjacent parts of the tank or escaping from the system (Fig. 1b).Stepwise merging, draining, and escaping of melt batches redistributes the gas in the system and modifies the gas pressure field in the tank, which may then trigger new merging and escaping events.In such a way a line of 'communication' exists between gas batches in the tank.
The continuous gas production and migration results in the development of a dynamic structure inside the tank which is continuously rearranged but remains statistically unchanged.The analysis of the data recorded during the experiment reveal power law size statistics of gas batches present in the system at a chosen timestep (Fig. 1c; Urtson & Soesoo 2007) and the volumes of escaping batches (Bons & van Milligen 2001).Due to discrete escape events, the total amount of the gas in the system fluctuates, producing a signal with a 1/f power spectrum (Bons & van Milligen 2001;Urtson & Soesoo 2007).The tank consists of two glass plates with 6.5 mm space between them.The tank is filled with fine-grained sand and water-sugar-yeast mixture.The processes in the tank were recorded by a digital photo camera and subsequently analysed using image analysis and statistical methods.(b) Merging of gas batches.Scale bar is 3 cm.(c) Cumulative distribution of gas batch sizes at a selected timestep during the ballistical transport mode of the system.The sizes were obtained by measuring the apparent area of batches on digital images.The data follow a power law N (> S) = S -D with a distribution exponent D = 0.5.(d) Dynamics of the gas escape in a 20-min interval and estimated batch size distribution exponents at some selected timesteps.Major gas escape events are marked by a sudden drop in the total batch area.The exponent of the gas batch size distribution, however, remains unaffected (from Urtson & Soesoo 2007).
The study of gas batch size distribution at selected timesteps shows that regardless of the extraction of large gas volumes, the batch size statistics remain unaffected.As most of the extracted gas volume is included in the few largest batches, removing of them does not affect the overall trend of the data and the fit to a power law (see Fig. 1c).As a result, the distribution exponent remains constant even when a large fraction of the gas is extracted (Fig. 1d).This proves the ability of the system to adapt to higher transport rates and discharge any additional gas without modifying its structure.The dynamic state of the system is characterized by only one number -the distribution exponent D.
As the experiment suggests, the transport of large gas volumes can be accomplished without the need for a steady-state open network of drainage paths.The liquid escape threshold (Vigneresse et al. 1996) can be overcome locally and temporarily; a high bulk fraction of the melt phase is not needed.The process of gas transport in the tank is a continuous rearrangement of local migration channels and accumulation pools, which disappear when the gas is extracted.
The analogue experiments described here contribute to the idea that melt segregation, transport, and subsequent accumulation is a highly dynamic stepwise and intermittent process.The observed phenomena such as a large number of interacting members (gas batches in the experiments), adaptation of the system to higher transport rates, and the emergence of power law statistics and 1/f fluctuations are characteristic of a self-organized critical system (Bak 1996;Bons & van Milligen 2001).

LEUCOSOME WIDTH STATISTICS IN MIGMATITES
The spatial distribution and width of leucosomes have been investigated in several studies (e.g.Tanner 1999;Marchildon & Brown 2003;Soesoo et al. 2004a;Brown 2007).The results are somewhat diverse -both power law (fractal) (Soesoo et al. 2004a;Bons et al. 2004) and non-fractal (Marchildon & Brown 2003) leucosome width statistics have been reported, as well as fractal spatial distribution of leucosomes (Tanner 1999).We provide here new leucosome data from several localities in addition to the data published previously in Soesoo et al. (2004a) and Bons et al. (2004).
Migmatites at the Montemor-o-Novo outcrop (38°38.12′N,8°12.54′E) in central Portugal belong to the Evora massif of the Ossa Morena zone.The highgrade metamorphism of the Seria Negra Group metasediments is probably of Variscan or Cadomian age (Pereira & Silva 2002).
The thickness of leucosomes and their spacing were measured along the axis of drill cores or along line traverses on outcrops using a tape ruler.The resolution of measurements was limited to 2 mm; leucosomes with thicknesses below this value were not counted as their number would be very likely underestimated and thus should not be included in the data.
When plotted on a bi-logarithmic graph with measured thickness on the -axis x and the number of leucosomes on the -axis, y the data follow a power law which is defined by a straight line with the distribution exponent being equal to its slope (Fig. 2).Different methodscumulative frequency and density distributions were used for plotting the power law data on the graph for comparison, as it helps to evaluate the quality of the data (Davy 1993).The cumulative distribution represents the number of leucosomes N whose thickness is greater than a given thickness : ( )  , where D is the distribution exponent.An alternative way of data presentation is the density distribution, which represents the number of objects belonging to an interval, divided by the interval length dh as ( ) .number of objects approaches 1 (see Davy 1993;Bonnet et al. 2001).Also, the deviation of the data from the power law trend at the small scale of the distribution due to insufficient sampling of the smallest objects (truncation effect) is exaggerated by the density distribution which helps to define the range of the validity of the power law and thus estimate the correct distribution exponents (Bonnet et al. 2001).The size of the bins where objects are distributed according to their size is critical as it determines the smoothness of the data (Davy 1993;Bonnet et al. 2001).In this study logarithmic binning was used for plotting the density distribution data.
In addition to leucosome width data, the fraction of leucosomes or the amount of melt during the late stage of melting was estimated by integration of leucosome thicknesses along the measuring line.This fraction rather describes the minimal amount of melt as a large fraction of melt may reside in the smallest leucosomes remaining under the resolution limit.The observed migmatite outcrops an drill core samples yielded minimal apparent melt fractions from 19.3 to 71% (Fig. 2).
As can be seen in Fig. 2, a power law is followed in most migmatites over at least two orders of magnitude on the length scale.The span of the distribution of leucosome thicknesses is limited due to several reasons.On the larger end of size range the cut-off occurs due to the limited extent of the outcrop or drill core compared to the whole magmatic system.As a result of this spatial restriction, larger objects present in the system at a lower frequency of occurrence may not be encountered on the measured traverse.At the lower end of the distribution the undersampling of the smallest objects due to limited resolution results in the deviation of the data from the power law trend, also known as truncation effect.Alternatively, it has been proposed that smaller objects can be undersampled due to their real absence in the section, as natural processes may set limits to the power law extent, thus the involved process may exclude the formation of smaller leucosomes (e.g.Bonnet et al. 2001).For the above reasons the data points suffering from the truncation effect were disregarded in the fitting procedure.On the data presented in this study the deviation from the power law trend occurs in some cases above the set resolution limit -already at 5 to 10 mm.Most probably the reason for this is the combination of the above-mentioned factors plus the effect of the sampling of apparent leucosome thicknesses.In several cases the placement of the measuring traverse was dictated by the orientation of the drill core or the erosional surface of the outcrop which can be at some angle with respect to leucosome layering.As apparent thicknesses are proportionally wider than real thicknesses, the resulting leucosome thickness distribution exponent remains unaffected.
It should be noted here that the distribution of leucosome widths and gas batch sizes in the analogue experiment are not directly comparable to the corresponding volume distributions as they represent oneand two-dimensional cuts through a three-dimensional system, respectively.Soesoo et al. (2004a) have derived the relationship between the distribution exponents of magmatic leucosome widths and respective melt batch volumes, where the expected range of the volume distribution exponents m between good accumulation and total dispersion 2 3 1 m < < (see also Bons et al. 2004) translates into 0 1 D < < for leucosome width distributions.
The comparison and combination of different plotting techniques allows validating the best fit of the data to a power law trend (Fig. 2).The scale invariant nature of the distribution of melt in former migmatite systems thus supports the concept of magma segregation, transport, and accumulation as a self-organized critical process.Despite the variability among different localities, no correlation can be drawn between the distribution exponent and host rock composition/metamorphic history (i.e.melting conditions and reactions) of the studied migmatites.As the mobility of the melt is a critical factor in such self-organized critical systems (Bons et al. 2004), the evolution and state of the system right before solidification may be determined by the tectonothermal situation at late stages of melting.

RELATIONSHIP BETWEEN TECTONIC FORCES AND PARTIAL MELT CHEMISTRY AS EXEMPLIFIED FROM THE NUMERICAL MODELLING
Many granitic intrusives show evidence of an origin from multiple sources (Leeman et al. 1990;Poli 1992;Elliott et al. 1997;Soesoo & Nicholls 1999) and complex magmatic processes (Chappell & White 1974;Clemens & Wall 1984;Soesoo 2000Soesoo , 2006)).The melts from different sources may completely mix to form a range of hybrid phases, or they may remain separate mingled phases: magmatic enclaves are the classical examples (Vernon et al. 1988;Poli & Tommasini 1991;Elburg 1996;Flinders & Clemens 1996;Perugini & Poli 2000, 2004;Perugini et al. 2003Perugini et al. , 2006)).Although it may be evident from the end product that mixing and mingling did take place, it is less clear how and where these processes occurred.Thus the composition of magmas depends, amongst other factors, on processes of partial melting and mixing of the resulting melts (Brown 2007).Although this has been widely recognized, it is difficult to quantify the relationship between tectonic forces and the composition of the melts generated under a concrete tectonic regime.
The key ideas for the numerical model which is able to investigate the relationship between tectonic forces and melt production and accumulation were already expressed a decade ago (Soesoo & Bons 1998, 1999).The aim of the developed numerical model is to investigate the effect of stepwise segregation and accumulation of melt batches during progressive melting of a source region, assuming that melt resides and moves solely in mobile hydrofractures or veins.In the model melt-filled veins are modelled as discrete spheres.It has been shown (Bons et al. 2004) that the melt pocket shape in the model will not change the results.Progressive melting is simulated in discrete time steps, typically 5000-10 000 steps.Melt batches are created in the space of a cube by increments of 0.01% of melting (all input parameters are changeable if needed).The melt movement, i.e. the mobility of hydrofractures, is simulated by moving all spheres each step over a distance that is randomly chosen between ± E (E is ranging between 0 and 1000 in the simulation; i.e. three orders of magnitude) in the x-, y-, and z-direction.In nature tectonic stress gradients may lead to propagation in any direction (depending on the stress field, one direction is preferred), while buoyancy forces would lead to upward propagation.A preferred upward propagation is not included in the model.This implies that the model applies to the region and scale where tectonic forces dominate over buoyancy driven melt transport.The value of E exemplifies the deformation rate.Merging the melt batches is performed by merging all pairs of overlapping or touching spheres where the new position of the sphere (melt batch) is taken as the volume-weighted average of the centres of the merged spheres.The extraction of melt batches is related to the batch volume: any sphere that exceeds a set threshold in volume (usually taken as accumulation of 2000 initial spheres or 1-1.5% of partial melt) is removed from the cube, simulating the escape of a melt batch from the source region.Statistics of the system and the extruded batches are recorded at every time step.More details of the mechanical part of the model can be obtained from Bons et al. (2004).
In nature the extruded spheres represent batches of melt that are emplaced at higher levels, either by dyking or by other ways of pluton/magma chamber formation.Each of these batches may have a different chemistry, due to local differences in source rock, differences in the melting time and regime, and accumulation history (Stephens 1992;Hobden et al. 1999;Soesoo 1999;Sobolev et al. 2000).The major element chemical composition of the melts generated at the different stages of partial melting can be calculated using the existing modelling packages (i.e.MELTS software developed by Ghiorso & Sack 1995).While knowing liquid/solid proportions during the melting steps at certain pressuretemperature (PT) conditions, the trace element distribution between solid and liquid phases can be calculated using conventional batch and fractional melting equations.Thus, the chemistry of each melt batch can be recorded for the major elements (using MELTS) and complementary trace element concentration (partition) at certain temperature.At the early stage of melting the forming melt batches have high concentrations of incompatible elements.Due to mixing with later batches (with less incompatible elements), the resulting melt batch will have a different chemistry which is depending on the value of E and the trace element distribution between the existing melt batches.When spheres merge, a volume-weighted average and volume-dependent average chemistry is calculated for each new batch, thus simulating the geochemical effect of mixing of different melt batches.
The results of the modelling show that smaller percentages of melt are needed to produce escaping melt batches (intrusives, dykes) at a high strain rate than at a low strain rate.The mobility factor E has a distinct effect on the character of the population of extruded spheres and their chemistry.With increasing E the average size of the extracted batches decreases (the average age also decreases), but the variability of ages increases.In geological terms this implies that a high deformation rate within the source would favour many, but smaller batches to leave the source, and a high variability in melt chemistry, reflecting the relatively short residence time within the source.Deformation distinctly facilitates the accumulation of melt batches that initially formed at different melting events during progressive melting and in different source sub-regions with variations in chemical composition.The early melts produced at a high strain rate are enriched in incompatible elements.The abundances of most incompatible elements may vary by two to three orders of magnitude between the early and late magma batches/dykes (Fig. 3).As melting proceeds, the extracted partial melts may show a large variation in chemical composition.At low strain rates the residence time of the initial melt batches is longer, enabling the melt batches to become equilibrated with residual solids and the chemical variation in extracted melts is small (Fig. 3).At a high strain rate (E = 1000) the trace element compositions of extracted melt batches match well with the theoretical incremental melting curve, while at a low strain rate (E = 2-10) the majority of melt batches are compositionally similar to those of theoretical batch melting.At a very high degree of partial melting (40-50% of source melting), the resulting melts are compositionally between these two end-member cases (Figs 3,4).Similar features can be noted if the concentrations of incompatible elements are plotted against a moderately compatible element as shown in Fig. 4.  Theoretical batch and incremental melting curves are for references.Note that at a high strain rate the concentrations of these elements in the melt follow the theoretical incremental melting curve, while at low strain rates the resulting melts show geochemical similarities with those for batch melting.
Apart from P, T, source composition, and partial melt percentage, differences in the deformation history and related differences in melt accumulation and extraction can also result in a large spread in geochemical characteristics of partial melts within the crust.This model shows that chemical variation within intrusive or extrusive complexes may not always be a result of different source rocks for partial melts or melting percentage alone but some variation can readily be caused by variation in tectonic stresses in the area.It can be shown that a relatively large change in the major and trace element compositions may be related to differences in strain rates.Figure 5 demonstrates results of a progressive water-saturated melting of a mixed metasedimentary source consisting of a 50 : 50 mixture of psammitic and pelitic rocks from the Proterozoic Halls Creek sedimentary sequence, Kimberley region, Australia.The source rocks were similar in major element composition, containing 71 and 73 wt% SiO 2 , respectively.At a low strain rate the resulting expelled partial melts have a very narrow range of SiO 2 contents (around 74 wt%), while at higher strain rates (10-100) the partial melts have a relatively large range of SiO 2 contents (72-77 wt%; Fig. 5).A similar relationship is valid for most of the modelled trace elements and several major elements.
One of the conclusions of this modelling is that the interaction between melt and source may not be simple and might not be described meaningfully by average values.Some melt batches have a long residence time, while others may quickly become entrained and have little time to interact with the restite or other melt.Melt chemistry may be variable from one extracted batch to another.In some cases the trace element, but also some of the major element chemistry can show large variations.The modelling results of accumulation efficiency also show that there is always some delay between melt formation and extraction.This implies that the volume of melts emplaced at a higher crustal level is always lower than the actual amount of melt generated in the source, especially when partial melting occurs during a tectonic quiescence period (low strain rate).
Since in our model melt escapes from the source before a percolation threshold or 'critical melt percentage' (Arzi 1978) has been reached, melt extraction starts at an earlier stage than in conventional models of granite This example uses a natural input data on the metamorphosed sequence of mixed psammitic and pelitic rocks.Both fractional and batch melting results (using MELTS software for major element and conventional equations for trace element calculations) for modelling the partial melting at 4 kb and subsequent melt accumulation and extrusion are displayed.Note the large range in SiO 2 contents at a high strain rate.The major element liquid compositions were calculated for every 3 degree centigrade, which generally corresponded to melt fraction of ca 1%.Trace elements were calculated according to the corresponding melt fractions.
petrogenesis.This will reduce the time for equilibration with restite or source rock, which in turn may cause disequilibrium partitioning of trace elements and isotopic ratios (Harris et al. 1995;Tommasini & Davies 1997;Perugini et al. 2006).For granitology it may indicate that granite geochemistry may not 'image their sources' (Chappell & White 1984) and attempts to model the composition of the source by equilibrium partitioning between melts and restite would also be doomed to fail (Sawyer 1991).However, this approach may lead to better explanations on the heterogeneous chemistry of some batholiths.The geochemical features of batholiths have often been explained by heterogeneous occurrences of restite phases or minerals or in situ differentiation (Michael 1984;Chappell 1996).Our analogue and numerical modelling results indicate the formation of large magmatic bodies by coalescence of distinct melt batches with somewhat different chemistry (see also Figs 3-5).Similar views have been emphasized by a number of researchers (e.g.Wareham et al. 1997;Pressley & Brown 1999;Slater et al. 2001;Perugini et al. 2003;Spiegelman & Kelemen 2003;Bons et al. 2004;Soesoo et al. 2004a;Urtson & Soesoo 2007).
Initially, the geological applications of the model were focussed on crustal melting and the generation of felsic melts (Soesoo & Nicholls 1999).However, this model can potentially be utilized in studies of melting of the shallow mantle, basaltic melt accumulation, and ascent into the crust.Moreover, this model can be readily developed for application to more complex geological histories such as a combination of both shallow mantle and crustal melting with subsequent mixing of melt batches within the framework of certain tectonic regimes, or even for hydrocarbon extraction (Bons & Soesoo 2003).The results derived from this model are likely to approach reality better than those from more conventional geochemical models (e.g.trace element modelling).

CONCLUSIONS
Different methods have been applied to the study of the behaviour of melts in partially molten systems, including numerical and analogue modelling techniques.The results contribute to the idea that the processes of melt production, its segregation from the solid matrix, and subsequent transport and accumulation are highly dynamic and stepwise processes and exhibit scale invariant patterns characteristic of self-organized critical systems -power law batch size statistics and 1/f volume fluctuations.The power law leucosome thickness statistics observed at several migmatite localities allow drawing parallels between the behaviour of analogue-experimental and natural migmatitic systems.To understand how magma composition evolves in the light of this new perspective, new models need to be developed that can take the stepwise character of the partial melting processes into account.
The numerical model presented in this paper enables qualitative and quantitative assessment of the effects of stress-induced melt migration and accumulation on the chemistry of partial melts.Results clearly show that at a lower strain rate it takes a longer time and more extensive melting of the source before the first melt batch leaves the source.A smaller percentage of melt is needed to produce extraction at a high strain rate.This in turn means that the residence time of small initial melt pockets is shorter at a high strain rate, which would result in liquid-solid disequilibrium and would influence the geochemistry of the extracted melts.
Geochemical variations within and between intrusive/ extrusive complexes are thus not only caused by different sources/melting percentages, but can also be due to other processes such as deformation.The present model can take into account the effects of different rates of melt extraction and accumulation on magma chemistry and may quantitatively contribute to the understanding of the relationship between tectonic environment and granite chemistry in contrast to static models applied previously.

Fig. 1 .
Fig. 1. (a)Setup of the experiment.The tank consists of two glass plates with 6.5 mm space between them.The tank is filled with fine-grained sand and water-sugar-yeast mixture.The processes in the tank were recorded by a digital photo camera and subsequently analysed using image analysis and statistical methods.(b) Merging of gas batches.Scale bar is 3 cm.(c) Cumulative distribution of gas batch sizes at a selected timestep during the ballistical transport mode of the system.The sizes were obtained by measuring the apparent area of batches on digital images.The data follow a power law N (> S) = S -D with a distribution exponent D = 0.5.(d) Dynamics of the gas escape in a 20-min interval and estimated batch size distribution exponents at some selected timesteps.Major gas escape events are marked by a sudden drop in the total batch area.The exponent of the gas batch size distribution, however, remains unaffected (fromUrtson & Soesoo 2007).
the exponent of the density distribution is increased by one compared to the exponent of cumulative distribution density D = cumulative 1 D + (Bonnet et al. 2001).The density distribution has the advantage of being free of the unfavourable curvature (censoring) effect of the cumulative distribution at the larger scale where the

Fig. 3 .
Fig.3.Melting step (melt%) vs. incompatible element concentration plot of extracted melt batches as a function of time during progressive melting.Theoretical batch and incremental melting curves are for references.The initial concentration (C 0 ) of the incompatible element (partition coefficient Kd = 0.01) in the source rock is 30 ppm.Note that at a high strain rate the concentration of the incompatible element in the melt follows the incremental melting curve.

Fig. 4 .
Fig.4.Incompatible element (Kd = 0.01) vs. moderately compatible (Kd = 2) element plot of melt batches derived during progressive melting.Theoretical batch and incremental melting curves are for references.Note that at a high strain rate the concentrations of these elements in the melt follow the theoretical incremental melting curve, while at low strain rates the resulting melts show geochemical similarities with those for batch melting.

Fig. 5 .
Fig. 5. SiO 2 contents versus temperature in the melting source (at the moment of extraction) and as a function of the strain rate.This example uses a natural input data on the metamorphosed sequence of mixed psammitic and pelitic rocks.Both fractional and batch melting results (using MELTS software for major element and conventional equations for trace element calculations) for modelling the partial melting at 4 kb and subsequent melt accumulation and extrusion are displayed.Note the large range in SiO 2 contents at a high strain rate.The major element liquid compositions were calculated for every 3 degree centigrade, which generally corresponded to melt fraction of ca 1%.Trace elements were calculated according to the corresponding melt fractions.