Double Power-Law in Leucosome Width Distribution: Implications for Recognizing Melt Movement in Migmatites

Leucosome sizes in migmatites has been shown to follow power-law distribution, which is indicative of self-organized criticality governing accumulation and transport of partial melts in anatexis. In our measurements of leucosome widths in the Olkiluoto migmatite in western Finland, we found double power-law distributions in addition to single power-law distributions. We divided the double power-law behaviors into those with positive and negative kinks separating the two power-law regimes. We propose that double power-law distributions in migmatites result from impediments on the bottom-up self-organized criticality that melt transport processes commonly display. These impediments can be caused by some leucosomes solidifying before others, in which case the leucosomes were effectively formed in two distinct events or phases. We also consider the emergence of double power-laws as a result of melt accumulation in and sudden melt removal from conduits of a specific size.


INTRODUCTION
Migmatites are high-temperature metamorphic rocks formed by anatexis (partial melting) in continental middle and lower crusts. As any rock will start to melt at sufficient pressure and temperature conditions, often aided by the influx of fluids, and the forming melts depend on the protolith (source rock) composition, migmatites can be very heterogeneous in chemical and mineralogical composition (Sawyer, 2008). The heterogeneity combined with their messy appearanceusually striped or flecked, and easily affected by any tectonic forcesmakes migmatites complicated to study. In addition, within the same or temporarily close tectonomagmatic cycle, several melting events may occur (e.g., Korhonen et al., 2010), and melt is transported from some parts of the anatectic area into other parts (D'Lemos et al., 1992;Brown, 2007). Effectively, migmatites are the crystallized remains of a melt factory, providing solidified snapshots of how partial melting in the crust looked at certain times. The complex melting and melt transport processes leading to that situation are, however, not easily identifiable from the field appearance.
Migmatites consist of paleosome, the metamorphic rock that did not melt, and neosome, the rock that was formed through anatexis (e.g., Sawyer, 2008). Furthermore, the neosome consists of melanosome, the rock from which the partial melt was extracted, and leucosome, an igneous rock that represents the crystallized partial melt. Leucosomes are useful as proxies for melt pockets in anatexis, because anatexis as a process is more difficult to observe than the crystallized migmatites it produces. However, it should be kept in mind that leucosomes only record the state of partial melts at the time of crystallization, which can differ from the period when the melt was active (Sawyer, 2001). Based on the proportional amounts of leucosome and paleosome, migmatites are commonly classified into leucosome-rich diatexites and paleosome-rich metatexites (Sawyer, 2008). This classification is descriptive and does not consider if the difference in melt proportions stems from higher grade of partial melting in diatexites than metatexites, or from melt transport processes that have depleted metatexites and enriched diatexites in melt.
Anatexis and transport of partial melts within an anatectic crust and outside it into granitic plutons is a significant process that contributes to the chemical differentiation of the continental crust and drives orogens to grow laterally (Jamieson et al., 2011;Sawyer et al., 2011). Ductile shear zones are known to act as melt transport channels in migmatite terrains (e.g., D'Lemos et al., 1992;Vigneresse, 1995;Brown and Solar, 1998;Johannes et al., 2003;Brown, 2007), although tectonic shear zones and partial melt transport are not always connected. Lee and others (2018) have shown that not all shear zones in an anatectic crust transport melt and, conversely, Bons and others (2009) have shown that some structures previously identified as shear bands and boudin necks are collapse structures formed by melt removal, not dependent on differential stress and thereby not tectonic shear zones. For an anatectic melt to move into shear zones and upwards in the crust, it must first be removed from the sites of partial melting and collected into a shear zone, necessitating melt transport processes that start from the grain-scale and continue through the outcrop-scale to the major ductile shear zone scale. Some authors see this melt transport in anatectic areas occurring through an interconnected network, where melt flows from small rivulets of melt into ever larger melt pathways, resulting in a network of connected melt-bearing veins of many sizes (Brown and Solar, 1998;Tanner, 1999;Weinberg, 1999;Sawyer, 2001;Brown, 2007;Vanderhaeghe, 2009). Others, like Bons and others (2009) propose that stable connections between melt accumulations are untenable, as local and transient transport of melt would cause the collapse of connectivity, and propose that melt accumulates in a stepwise manner. According to this model, melt-bearing veins or pockets only make temporary connections to each other when melt moves in batches from one site to another, transport pathways closing after a melt batch has moved. What most of the models of anatectic melt transport have in common is that they view it as a bottom-up process that exhibits self-organization.
Self-organized criticality (SOC) is a process where statistical order is generated from many individually unpredictable events in a system with both spatial and temporal degrees of freedom (Bak et al., 1988). Systems displaying SOC behavior follow powerlaw distributions, and therefore a power-law distribution can indicate SOCalthough it should be noted that other processes can also produce power-law distributions (Corral and González, 2019 and references therein). A power-law distribution is often visualized as in Figures 1A,B: either a negatively sloping straight line or a downward curving line on a log-log plot of the cumulative number of events as a function of event size, here shown for leucosome widths. In addition to this single line distribution, other forms of power-law distribution also exist. A double (or broken) power-law distribution consists of two power-law distributions with different slopes, separated by a kink at a specific event size d as shown in Figures 1C,D. In this work we refer to the general shapes of power-law distributions in Figure 1 as single untruncated ( Figure 1A) or truncated ( Figure 1B) power-law distributions, and as double power-law distributions with a positive kink ( Figure 1C) or a negative kink ( Figure 1D). In Earth sciences, double power-law distributions have been shown to describe phenomena such as the severity of tornadoes in the United States (Pinto et al., 2014), the size of FIGURE 1 | Simplified sketch of the different leucosome width distribution types. (A) In migmatites showing this leucosome distribution pattern, leucosome widths follow a single untruncated power-law distribution. The value of m indicates the slope of the line: if m is above 1, most of the solidified melt resides in small leucosomes, but if m is close to 2/3, the largest leucosomes or dikes host the majority of the solidified melt. (B) The leucosome widths follow a single truncated power-law distribution, resulting in the curved shape of the function. Truncation can be an artifact from measuring, in which case it indicates that the measured section was not long enough to contain the largest leucosomes. Alternatively, truncation can result from the power-law distribution being limited, in which case it indicates that the largest measured leucosomes are nearly in the same order of magnitude as the largest possible leucosomes that follow this power-law distribution. (C) The smallest leucosomes follow a different power-law distribution from the largest leucosomes. The change from the first trend to the second occurs at around width d. The leucosomes close to width d are relatively more abundant than in (A). (D) The smallest and largest leucosomes follow different trends, but the situation is the opposite of (B) most of the melt is found in the very smallest and very largest leucosomes, while leucosomes sized around width d are underrepresented.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 591871 forest fires in Portugal (Pinto et al., 2014), the daily flow above the annual flow in some hydrological systems (Segura and Pitlick, 2010;Segura et al., 2013), and the seismic moment of global earthquakes (Corral and González, 2019). Many studies concerning power-law distribution of leucosomes in migmatites have been presented (see chapter two for a literature review), but the focus has been on identifying single power-law trends and establishing that partial melting displays SOC. In this study, we measured leucosome width distributions and identified double power-law trends in them to show that the SOC behavior of partial melting can, in some cases, result in leucosome sizes that do not follow a single power-law distribution. We chose Olkiluoto Island in western Finland as our study area, because in addition to hosting a suitable migmatite, the metamorphic history of the area is well constrained, and the location of shear zones is known. We compare the leucosome size distributions in different parts of this spatially limited migmatite area to one another and to previous geological research on the anatectic history of the area in order to estimate where double power-law distributions stem from. We model numerically two sets of partial melts that display SOC independently of one another, to quantify how leucosome size distributions would appear on a single outcrop where two partial melting events occurred. We also consider other mechanisms that may cause double power-law distributions in leucosome widths. Our findings can be applied to identify multiple sets of leucosomes on other migmatite areas where leucosome sizes do not follow single power-law distributions and aid in constructing the anatectic history of such migmatites.
Power-law distribution of leucosome sizes in migmatites can be seen as the probability of encountering a leucosome of a certain size or larger. The smallest leucosomes are ubiquitous, but the larger the leucosome, the less likely it is to have formed. This means that the probability of any given leucosome in a migmatite being small rather than large is very high. This is expressed as: where N > x is the number of observed leucosomes that are larger than size x (width in mm), k is a constant defining the total amount of leucosomes larger than or equal to unit size, and m the distribution exponent that gives the behavior of the distribution over different orders of magnitude. The exponent m is an empirical constant whereas the unit size can be chosen arbitrarily, taken to be 1 mm in this work. When the function is plotted on a log-log graph, the result is a straight line with the slope defined by m and the value at unit size given by k, as shown in Figure 1A.
Power-law distributions in nature are always valid only for a limited range of values for the observed feature. Truncated power-law distributions are used to describe situations where the largest events either do not follow the same power-law distributions as the smaller events or where large events have not been observed due to measurement limitations (Corral and González, 2019). In some cases, two different power-law distributions can be identified in one data set so that events smaller than a cutoff size d follow a truncated power-law, and events larger than d follow another power-law distribution (Corral and González, 2019). These are called double power-law distributions.
The exponent m in the equation describes if leucosomes in a migmatiteand by extension, melt-bearing pockets or veins during the anatexisare predominantly small or large. When m is low, the line extends far toward large leucosome sizes, meaning that there are proportionally many large leucosomes. The higher value m takes, the steeper the line and the more numerous the small leucosomes are. In other words, most of the melt was concentrated in large conduits if m is low, and if m is high, most of the melt was situated in many small pockets . Typical values for m in leucosome distributions range from about 2/3 to 1.9 (Urtson and Soesoo, 2009;Yakymchuk et al., 2013;Soesoo and Bons, 2015).
Many leucosome size distribution studies have noted that the power-law functions that describe the leucosome size distribution for most of the size range do not apply for the smallest leucosomes (Bonamici and Duebendorfer, 2010;Yakymchuk et al., 2013). In many cases, these discontinuities in leucosome distribution diagrams have been attributed to measurement limitations (Bonamici and Duebendorfer, 2010) or size range limits for the power-law behavior (Marchildon and Brown, 2003;Yakymchuk et al., 2013). Marchildon and Brown (2003) suggest that the original anisotropy in partially melted metasediments can cause this limitation, but it has also been found in Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 591871 metaigneous rocks that lack the strong anisotropy of metasedimentary rocks (Yakymchuk et al., 2013;Saukko, 2016). When one measures line transects of leucosomes, measurement limitation will affect both the smallest and largest leucosomes. These measurement limitation effects are visible in the distribution graphs as deviations from the power-law trends. When logging leucosomes of the minimum measurement resolution width, it is difficult to determine whether a given leucosome is large enough to be included in the measurement set, or small enough to be left out, leading to possible misclassification at the size range. Measurement limits can also affect the largest measured leucosomes, since the probability of a very large (i.e., uncommon) leucosome occurring on a measured section is effectively random: in some sections, there will be large outliers while they are absent on others (Bonamici and Duebendorfer, 2010). Thus, measurement limitations affect the number of leucosomes close to the minimum measurement resolution width and the largest leucosomes, but significant misclassification is unlikely to occur between these sizes. Any error in the number of large leucosomes affects the whole distribution due to it being cumulative, but the relative size of the error is small for the more numerous, smaller leucosomes.
We propose that kinks in leucosome width distribution diagrams can signify underlying double power-law distributions: leucosomes smaller than the kink size follow a power-law distribution that is different from the power-law that governs the size distribution of the larger leucosomes. Thus, these kinks are not necessarily caused by measurement limitations or truncation. Multiple power-law regimes within size distribution diagrams have previously been identified by Soesoo and Bons (2015) and Saukko (2016).
As mentioned in the introduction, many models of melt transport are based on self-organized criticality and predict power-law behaviors. Therefore, field observations of powerlaw behaviors displaying features such as kinks or very high or very low exponents, can provide valuable insights into the underlying processes of melt accumulation and transport. In a SOC state (e.g., Bak et al., 1988), a system adjusts itself to accommodate transport, whereas in a classical Darcian system, the transport adjusts itself to the system . For magma transport, this means that the rock matrix is modified to allow transport to happen at the rate determined by the melting rate, and the amount of magma residing in melt pockets and veins. A typical feature of a SOC system is that there is a strongly dynamic balance between input and output and that any small perturbation canbut not necessarily mustlead to a large chain reaction or an avalanche. In magma transport, this "avalanche" would be a large-scale merging of melt batches and emplacement of large melt volumes.
According to the stepwise melt accumulation model proposed by Bons and others Bons et al., 2009), Bons and van Milligen (2001), Soesoo and others (2004), and Urtson and Soesoo (2007), magma batches join each other in a stepwise manner, continuously merging into ever-bigger batches and migrating farther away from their point of origin. Eventually, batches large enough to escape the source are formed, and magma can ascend through the crust as dikes to feed larger igneous bodies. Each time two or more veins are connected, transport can occur and lead to the merging of the magma in these veins, but also to the collapse of a part of the veins involved and, thus, the destruction of local connectivity (Bons et al., 2009). Contrary to classical models of magma percolating through a network of everbigger veins (e.g., Weinberg, 1999), veins are rarely connected, as their connection is only transient. This model is supported by analogue experiments (Bons and van Milligen, 2001;Urtson and Soesoo, 2007) and by numerical modeling . The observation in outcrops that leucosomes tend to be parallel, with generally few visible second order melt pathways that crosscut the smaller leucosomes has also been suggested to be evidence for transient connections (Soesoo and Bons, 2015).
In areas with high melt production, melt can accumulate in and escape as batches multiple times, like air behaves in an accordion (Yakymchuk et al., 2013). Not all melt is necessarily lost when a batch empties-partial crystallization and the retention capacity of the rock affect the amount of melt that will stay trapped (Yakymchuk et al., 2013). Due to this, smallsized leucosomes are not necessarily indicative of small degree of melting: instead, they can indicate that melt has been lost from the migmatite. Many researchers (e.g., Vigneresse et al., 1996;Yakymchuk et al., 2013) have proposed that there is a melt extraction threshold, meaning a maximum size range of melt accumulation, at which it will empty instead of growing more. The volume or spatial dimensions of melt at the threshold are not known, and it is also unclear if similar melt accumulation thresholds occur at larger scales of melt transport, i.e., if they also show scale invariance.
Melt extraction efficiency can be roughly estimated from leucosome size distributions . The exponent m in Eq. 1 that describes how melt is situated in a migmatite can show if a pressure minimum appearing in the proximity can be effectively filled by melt extraction from the system. When m is low, a larger portion of melt resides in large melt-bearing veins rather than in small pockets. Melt transport over larger distances is restricted when the volume of melt in a system is below about 20% (Vigneresse et al., 1996) and an interconnectedness of melt is needed for melt migration (Vanderhaeghe, 2009), indicating that small and separate melt pockets may not empty as effectively as more voluminous collections of melts. Thus, melt can be removed more efficiently from a partially melted rock if it is situated in a single vein, whereas it would require more disruptions in the surrounding rock to empty as much melt from numerous small melt-bearing veins or pockets that are separate from one another.

GEOLOGY OF THE STUDY AREA
This study was conducted on the six by 2.5 km large Olkiluoto Island on the west coast of Finland, where the bedrock is particularly well studied (e.g., Tuisku & Kärki, 2010;Aaltonen et al., 2016) due to plans for using the area as an underground nuclear waste repository. The bedrock in southern and western Finland was formed in the Precambrian Svecofennian orogeny, at Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 591871 around 1.9-1.8 Ga (e.g., (Lahtinen et al., 2005;Korja et al., 2006;Nironen et al., 2017). During the orogeny, igneous rocks in the region formed in two distinct compressional phases, at 1.89-1.87 and 1.84-1.80 Ga (Väisänen et al., 2002;Lahtinen et al., 2005), and an extensional phase between them, at 1.87-1.84 Ga (Lahtinen et al., 2005;Korja et al., 2006). The bedrock of Olkiluoto predominantly consists of metatexitic migmatite in the western part and diatexitic migmatite in the eastern part ( Figure 2; Aaltonen et al., 2016). The early stages of the Svecofennian orogeny are recorded in the parent material of the migmatites: the metapelitic, meta-arenitic and intermediate metavolcanic rocks were deposited at 1.90-1.88 Ga (Aaltonen et al., 2016). The earliest ductile deformation stage discovered on Olkiluoto, D1, also occurred at this time (Aaltonen et al., 2016). A metamorphic peak at about 6 kbar that occurred at some time between 1.88 and 1.86 Ga was largely overprinted by later events (Tuisku and Kärki, 2010;Aaltonen et al., 2016). Starting at 1.86 Ga, during the regional extensional event, a series of tonalitic, granodioritic, and granitic rocks formed in Olkiluoto (Aaltonen et al., 2016). Deformation stage D2 and the second metamorphic peak, which erased much of the traces of the earlier one, occurred during the later regional compressional phase at 1.84-1.83 Ga, at 660-700°C and 3.7-4.2 kbar (Tuisku and Kärki, 2010;Aaltonen et al., 2016). D3 followed the metamorphic peak at around 1.82 Ga, and the final deformation stage D4 at 1.81-1.79 Ga. Production of granitic melts was at its most voluminous during the metamorphic peak at ca 1.84 Ga, but also occurred during the later ductile deformation stages D3 and D4, so that the last zircons in granitic veins did not crystallize until 1.80 Ga (Mänttäri et al., 2010;Aaltonen et al., 2016). Migmatization also occurred earlier, coevally with the granitoid magmatism at 1.86 Ga, as evidenced by zircon ages in leucosome pegmatoids (Aaltonen et al., 2016) The error margins of the U-Pb data overlap (Aaltonen et al., 2016), making it unclear if partial melts were present during the whole period from 1.86 to 1.80 Ga, or if there were multiple separate anatectic events. Further uncertainties concern the earlier metamorphic peak at about 6 kbar: establishing its timing and whether migmatization occurred during it has been obfuscated by the later metamorphic peak.  Koistinen et al., 2001). Main picture: simplified bedrock map of Olkiluoto, after Aaltonen et al., 2016. In general, the Olkiluoto migmatite is more diatexitic to the east and metatexitic to the west, although the transition is more gradual than indicated in the map. Strain was concentrated into the ductile shear zones that divide the island into several less-intensely deformed structural blocks. Shear zones oriented E to W are the oldest, and the younger the shear zone, the more NNE-SSW orientation it has. The measured drillcore sections are projected to the corresponding structural units on the surface according to Hartley et al. (2017). Labels A-D correspond to the sections in The paleosomes of the Olkiluoto migmatite were originally turbidite-type sediments, volcanic deposits, and felsic intrusive rocks (Aaltonen et al., 2016). These rocks now appear as different kinds of gneisses migmatized to varying degrees. The metapelitic gneisses on Olkiluoto typically consist of quartz, plagioclase, K-feldspar and biotite as major phases, and accessory cordierite and sillimanite (Aaltonen et al., 2016). The metapsammitic gneisses are simpler in composition, consisting of quartz, plagioclase, and K-feldspar as major phases, and biotite as a minor phase (Aaltonen et al., 2016).
Petrographically, leucosomes on Olkiluoto can be divided into tonalitic leucosomes consisting mainly of quartz and plagioclase with small amounts of potassium feldspar, and granitic leucosomes that consist of potassium feldspar, quartz, and plagioclase (Aaltonen et al., 2010). As the granitic leucosomes cross-cut the tonalitic leucosomes on some outcrops, they have been considered to be younger (Paulamäki and Koistinen, 1991;Kärki, 2015). On many other locations, no cross-cutting relations have been observed, and no geochronological studies show the absolute age of the different leucosomes in one migmatite. It is therefore unclear if the cross-cutting implies a significant age difference between the leucosomes, or if it simply shows that the tonalitic leucosomes were products of early crystallization, whereas the granitic leucosomes represent liquids differentiated from them.
The long tectonic history is visible in the bedrock of Olkiluoto Island also in the form of several shear zones that divide the island into less-deformed structural blocks (Figure 2; Aaltonen et al., 2016). The shear zones, which are well described by Aaltonen and others (2016), can be tens of meters wide but vary in the outcrop scale from narrow shear bands to several meters wide zones. The oldest shear zone systems were formed during D2 but extensively overprinted during D3. They are dextral, strike E-W and dip moderately to steeply to the S (Aaltonen et al., 2016). A younger, purely D3 stage shear zone strikes NE-SW and dips moderately to the SE, and the youngest shear zones from D4 strike NNE-SSE and dip moderately to steeply to the ESE (Aaltonen et al., 2016).

MATERIALS AND METHODS
We made field observations of the Olkiluoto migmatite on outcrops and examined a drill core. On all observed sites, leucosomes that display variation in color and texture were present. The leucosomes can broadly be categorized into two groups: gray fineand even-grained leucosomes (grain size <2 mm), and white leucosomes that are generally coarser and more uneven in grain-size than the gray ones. The gray leucosomes correspond to the tonalitic leucosomes described by Aaltonen et al. (2010), and the lighter leucosomes to the granitic leucosomes described in the same work. Examples of each group in the same sections are shown in Figure 3. Some leucosomes display intermediate colors and grain-sizes, obstructing definite categorization. We did not observe unambiguous cross-cutting relations that would reveal the relative ages of the two leucosome sets, but earlier studies suggest that the white, coarser-grained leucosomes solidified later than the gray, even-grained leucosomes (Kärki, 2015).

Measurements
We measured leucosome width distribution in stromatic migmatites on one surface outcrop and in three different sections of a drill core. Measurements were made by placing a measuring tape perpendicular to the rock's schistosity and logging the leucosome thicknesses along the tape. The smallest measured leucosomes were 2 mm wide, anything smaller than that was left unmeasured. Similar methods for leucosome measurements have been used by previous researchers (e.g., Soesoo et al., 2004;Urtson and Soesoo, 2009). In the case of the drill core, which cuts the general schistosity of the rock at a 45°angle, corrections from apparent thickness to actual thickness were made with the sine rule. With the outcrop measurements, corrections were not necessary, as the schistosity is perpendicular to the measured horizontal surface. Because the length of the measured sections can be critical in whether power-law distributions are seen in the measurement results (Yakymchuk et al., 2013), we chose section lengths of approximately 5 m or longer. The measured outcrop consists of stromatic migmatite and is located by a D2 shear zone that was overprinted during the D3 stage, at around 1.83-1.82 Ga (Aaltonen et al., 2016). A total of 693 leucosomes were measured in the 10.1 m long section A in the outcrop. All leucosomes on the outcrop were nearly parallel, so no semi-perpendicular leucosomes were neglected by the measurement method.
The 870 m long drill core OL-KR29 penetrates many of the different structural units of the island (Hartley et al., 2017). The amount of leucosome in it varies: some parts are diatexitic whereas others are more metatexitic. Our measurements are from stromatic parts in the drill core. Drill core section B is from the structural block in the center of the island, not associated with any of the known local shear zones. The measured section is 9.4 m long, and we observed 162 leucosomes in it. Measured section C is from the middle of a D3 shear zone, which is a younger shear zone than the one represented on section A (Aaltonen et al., 2016). We observed 114 leucosomes in the 6.0 m long section C. Drill core section D, which is 4.9 m long with 152 leucosome observations in it, is from the structural block in the northern part of the island. This section is not associated with shear zones.

Fitting of Power-law Distribution on Measured Data
We fitted both truncated and untruncated power laws to our measured cumulative leucosome distribution. To account for measurement limitations, power-law distributions were not fitted to the smallest leucosomes. In cases where the largest leucosomes clearly deviated from the overall trend line, they were also disregarded. We determined the value of m using the maximum likelihood method, which gives the m value of the distribution that would be most likely to produce the measured data. This method has been shown to give good estimates of m for power-law distributions (White et al., 2008).
The log-likelihood for an untruncated power law is given by where a and b are the lower and upper cutoff size and g the geometric mean of the leucosome width data (Deluca and Corral, 2013). The geometric mean is calculated as where N is the total number of measured widths and x i the measured widths. This means that the fits are made to the measured widths, yielding the slope of the cumulative distribution m. In the case of an untruncated distribution, the maximum of the likelihood function can be calculated directly, yielding For the untruncated power laws, we chose a value for the lower cutoff, which in combination with the measured leucosome sizes immediately gives a value for m. We then chose the value for k so that the fitted distribution matched the measured distribution for the lower cutoff size. To use maximum likelihood for fitting truncated power laws, first both the lower, and upper cutoff must be chosen and then the likelihood function must be maximized numerically (Deluca and Corral, 2013). We used the GRG nonlinear solver in microsoft excel for the numerical maximization and obtained the m value. We chose the value for k and shifted the fitted distribution so that it matched the measured distribution for both the lower and upper cutoff size (Corral and González, 2019). To ensure that the fits follow power-laws instead of log-normal distributions, we also tried fitting log-normal distributions to sections B, C and D using the microsoft excel function LOGNORM. These fits are shown in the Supplementary Material. Finally, we performed the Kolmogorov-Smirnov test to determine the best fit for each measured section.

Distribution Modeling
To better understand the double power-law behavior, we modeled how leucosome width distribution would appear on a hypothetical section if there were two sets of leucosomes, each following their own power-law distributions. Measurements from such section would simply represent the sum of leucosomes originating from both events. We let one synthetic set of leucosomes follow an untruncated power-law distribution with a low value for the exponent, m 2/3, and another an untruncated power-law distribution with a high exponent, m 1.9. We calculated the expected number of leucosomes for some sizes for each of the two distributions, summed them up, and plotted the new leucosome size distributions. We used microsoft excel to fit power-laws to these synthetic width distributions to approximate the resulting double power-laws.

Measurements
Leucosome width distribution results, shown in Figure 4, reveal that distribution types vary within the island of Olkiluoto. Two of Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 591871 our sections follow kink-free single power-law distributions and two display kinks indicative of double power-law distribution. Data from the outcrop section A show a negative kink at width d 24 mm, whereas the drill core section B displays a positive kink at width d 20 mm. Leucosomes measured in the drill core section C, from the middle of a D3 shear zone, follow a single truncated distribution. Drill-core section D from outside of shear zones in the northern structural block displays a single untruncated distribution, although the shorter section length may affect the representativeness of these results and can explain the unexpected distribution of large leucosomes on the section.

Modeling
Our numerical modeling of untruncated power-law width distributions in synthetic leucosome sets shows that distinguishing between two different distributions is possible in some cases, whereas in others the sum simply looks like a single distribution. Figure 5 shows a double power-law distribution with a negative kink for synthetic leucosome widths, akin to the measured outcrop section of Figure 4A Additional graphs illustrating the modeling can be found in the Supplementary Material. A negative kink emerges when the number of mid-sized leucosomes originating from each of the two original distribution sets are roughly equally numerous, and the exponents for the two distributions have very different values. The sum of the two sets then follows the original distribution with higher m at small leucosome widths, and the original distribution with lower m at large leucosome widths. If the leucosomes of one set are much more numerous than those of the other, or when both sets have very similar distribution exponents, the sum of the two individual distributions looks like a single distribution. Consequently, identification of more than two sets of leucosomes is difficult: for more than two kinks to be visible, all the sets should contribute with a large enough fraction of all leucosomes, and their power-law distributions should all have significantly different values for m and k. Because the leucosome set containing more numerous leucosomes in a given size range always dominates over the other set in that range when the sets are summed up, a positive kink in the width distribution graphs cannot result from the sum of two different untruncated distributions. A truncated distribution with low m and a lower cutoff at kink size d together with an untruncated distribution with high m and an upper cutoff at kink size d would form a positive kink, but to our knowledge there are no mechanisms in anatectic systems that could plausibly explain the emergence of the upper cutoff size for the distribution with the higher m.

DISCUSSION
The Olkiluoto Island is a spatially limited migmatite area, where the main aspects controlling anatexis, such as temperature, lithostatic pressure and regional stress direction, can be assumed to have been uniform. The discovery of multiple distribution types and double power law distributions in leucosome widths therefore reflects not the regional scale constraints but the inherent heterogeneity of migmatites.
Assuming that rocks undergoing anatexis are a self-organized system spanning over all the orders of magnitude that are observed in outcrop-scale measurements, leucosome widths should display single power-law distributions (e.g., Yakymchuk et al., 2013). While this often is the case both in our measurements and in previous research Soesoo et al., 2004;Urtson and Soesoo, 2007;Bons et al., 2009;Bonamici and Duebendorfer, 2010;Hall and Kisters, 2012;Yakymchuk et al., 2013;Soesoo and Bons, 2015), the occasional double power-law distributions show that sometimes different size ranges follow different distributions.
In the migmatites where leucosome widths follow double power-law distributions, the measured leucosomes are not parts of a single self-organized system, but could stem from multiple self-organized sets that were formed at separate times so that one forms on top of the other. Alternatively, the two systems may have been active simultaneously but only interacted intermittently. We present two different conceptual models for how these situations can arise in migmatites.

Multiple Melt Generations
Distribution type differences in the Olkiluoto migmatite mainly depend on the inherent heterogeneity of migmatites: because melts tend to drain from their partially molten sources once the melt volume is sufficiently high, some parts of crystallizing migmatites will represent areas where melts were drained from, whereas others represent sites where melts accumulated (Bons et al., 2008). As ductile shear zones are known to act as melt channels (e.g., D'Lemos et al., 1992;Vigneresse, 1995;Brown and Solar, 1998;Johannes et al., 2003;Brown, 2007;Bons et al., 2008), they can represent the migmatite parts where melt once accumulated in, and possibly later drained from into more extensive granite bodies higher up in the crust (Bons et al., 2008). Conversely, migmatite parts outside of shear zones are where the melts were formed in and where they drained from into the shear zones.
Felsic magmas in the Olkiluoto migmatite crystallized between about 1.86 and 1.81 Ga (Mänttäri et al., 2010;Aaltonen et al., 2016), so age differences between leucosomes can be large enough for some to have been crystallized while others were still forming (Kärki, 2015). The age difference between the melting events or melt pulses need not be as large as the extremes of the age range: fractional crystallization occurs simultaneously with transport of partial melts in migmatite terrains and partially crystallized leucosomes are significantly more competent than meltbearing conduits (Sawyer, 1987), suggesting that leucosomes undergoing fractional crystallization become less active as melt conduits. Additionally, during a long anatectic event, melt likely segregates from its source several times, so that leucosomes crosscutting each other can result from a later melt loss event creating new pathways that cross-cut the pathways formed in an earlier melt loss event (Diener and Fagereng, 2014). We propose that the double power-law distributions seen in two of the measured sections can be explained by one early set of melt-bearing veins and pockets that behaved as a self-organized system Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 591871 before crystallization (either complete or fractional) that rendered it inactive, and one later set of melts that formed its own selforganized system on top of the earlier, inactive one. Measurements from a D3 and a D2 shear zone plot as different distribution types. In the younger shear zone, we found a single set of leucosomes that fits a single truncated power-law distribution with m 0.76. In contrast, measurements from the older shear zone show a negative kink that separates the trends into two power-law distributions. The shear zone may have acted as melt channel when it was formed during D2 and expelled much of the melt in it into a higher crustal level. Thus, only the parts of melt-bearing veins that crystallized early or melts that were otherwise difficult to move remained at the site. These remains are now visible as the many thin leucosomes that result in the observed truncated power-law distribution with exponent m 1.78. Later, more melts accumulated in the shear zone, forming a collection of melts that solidified before it could be transported away. This melt pulse is seen as the set of leucosomes that fits the untruncated power-law distribution with m 0.95. Our numerical modeling of synthetic leucosomes fits with this model, showing that negative kinks can arise when two leucosome sets solidified independently of one another.
The age difference between the two leucosome sets in the D2 shear zone cannot be estimated from the distribution results or the mathematical modeling. The later melt pulse may have entered the shear zone during D2 but after the solidification of the first set of leucosomes, or later, during D3 when the shear zone was reactivated. In the latter case, the later melt pulse would be roughly coeval with the leucosomes in the D3 shear zone, which also display a comparable power-law exponent m between 2/3 and 1, the maximum and minimum melt transport efficiency values as established by Soesoo et al. (2004). If the leucosome sets with similar power-law exponents indeed are coeval, our results are in accordance with earlier assessments that less melt was Outcrop section from a D2 shear zone follows a double power-law with a negative kink, showing an underrepresentation of leucosomes around the kink size d (24 mm) in comparison to both smaller and larger leucosomes. (B) Drill core section 1 follows a double power-law with a positive kink (d 20 mm), corresponding to an overrepresentation of leucosomes of that size. This section is not from a shear zone. (C) A kink-free truncated power-law distribution fits the measurements from drill core section 2. This section is from the central parts of a D3 shear zone. (D) Drill core section 3 appears to follow an untruncated power-law distribution. This section is not close to a shear zone. The Supplementary Material shows the other distributions considered for the measurement data.
Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 591871 active in the anatexis during D3 than D2 (Aaltonen et al., 2016): during D3, there was enough melt to accumulate in the shear zones but not enough for it to move upwards in the crust, whereas melt removal had occurred earlier in the D2 shear zone. Similarly to the D2 shear zone, a new melting event in a migmatite where the earlier leucosomes were already at least mostly solidified may have caused the positive kink in the section that was measured outside of any shear zones. In this case, the first set of leucosomes would have formed in situ in the migmatite and the anatectic system would have formed a melt conduit size distribution following a single power-law. If melting proceeded long enough, the melt pockets would have grown and many of the largest melts could have escaped the rock after the system reached the melt escape threshold (Vigneresse et al., 1996;Rosenberg, 2001;Rosenberg and Handy, 2005), resulting in a truncated power-law distribution with low exponent m. When bedrock partially melts, the first melts are formed in the most melt-fertile locations, namely three-grain junctions commonly containing feldspar, quartz, and a hydrate phase (Brown, 2008). Such triple junctions frequently occur in metapelitic rocks, but because the most fertile locations produce melt first (Sawyer, 2008), fewer such points would be left in the metapelite for a second melt generation event. However, at the 3.7-4.2 kbar pressure and 660-700°C temperature conditions of Olkiluoto (Tuisku and Kärki, 2010), a neosome consisting of plagioclase and quartz in the leucosome and biotite in the melanosome could have partially melted if water was present (e.g., Weinberg and Hasalová, 2015;Sola et al., 2017). Thus, a second melting event in the Olkiluoto migmatite could have favored generating melts inside the previously formed neosomes rather than in the paleosome, where melt-fertile sites did not readily occur. This would have led into a situation where many of the older leucosomes were effectively re-melted, and due to local heterogeneities in the neosome composition, some could have been re-melted to a larger extent than others. Inheriting an earlier truncated leucosome distribution would have affected how the younger melts were distributed in the early stages of the second anatexis, possibly exaggerating the bend of the truncated distribution into a prominent kink. Given that melt production was lower during the later deformation stages, the second melt generation may well have stalled before it could reach self-organized criticality.
Field observations also support that there are two groups of leucosomes on the outcrop and the drill core sections (Figure 3). On the outcrop section, the majority of the observed leucosomes larger than the kink size 24 mm are white and coarse-grained, even though the vast majority of all measured leucosomes are gray and fine-grained. The gray leucosomes are generally small in all the drill core sections as well, even though the white leucosomes are more abundant than the gray on each drill core section, thus contrasting with the outcrop measurements. Earlier research also suggests that granitic leucosomes are generally more abundant than tonalitic leucosomes on Olkiluoto (Aaltonen et al., 2010). The differences between the two leucosome sets indicate differences in protolith compositions or in the crystallization processes between large and small leucosomes. This can stem from, for example, the differing nucleation rates in early and late melts (Marsh, 2015), melt fractionation caused by partial crystallization, or differing mineralogical, and chemical compositions of protoliths: the protolith of a partial melt formed at a later stage is inevitably somewhat different from the protolith of the earliest melts (Sawyer, 2008), even though the compositions of a normal metapelite and a metapelite modified by low-grade anatexis are not necessarily very large. If the first melts were formed in a metapelite but later melts in the tonalitic leucosomes, as we suggested above, the mineralogical difference of the protoliths is even larger.
Distribution measurements alone cannot show how the origins of the two leucosome sets on Olkiluoto differed from each other. Closer observations of the different types of leucosomes, especially detailed inspection of their crosscutting relations, could help ascertain if the leucosomes originate from separate anatectic events within the same migmatite, or if they represent different melt segregation pulses within the same anatectic event. Tuisku and Kärki (2010) found traces of an earlier metamorphic peak in a study of Olkiluoto migmatites, but it is not known if anatexis was associated with the metamorphic peak. Isotope-based mineral age determinations of zircon might help determine FIGURE 5 | Modeling results of leucosome width distribution when two untruncated power-law trends interfere with each other, comparable to the measured outcrop section ( Figure 4A). This double power-law distribution with a negative kink can form when two sets of leucosomes following different width distributions interfere. For the kink to be visible, the two distributions must have substantially different exponents m, and the total amount of leucosomes on the trend with the lower m value must be considerably lower than on the trend with the higher m value.

Melt Accumulation to and Melt Loss From "Kink-Sized" Melt Conduits
Another possible explanation of the kinks in leucosome distribution is that they reflect the processes of melt accumulation and loss in the migmatites. Yakymchuck and others (2013) suggested that melt accumulation and extraction is a cyclical process, where melt accumulates on the site of partial melting, escapes in a melt loss Frontiers in Earth Science | www.frontiersin.org November 2020 | Volume 8 | Article 591871 event once the melt extraction threshold is achieved, and starts accumulating again after the event. The size range of the threshold depends on the amount of melt and the strength of the melting crust (Yakymchuk et al., 2013). Solidified migmatites show an active system frozen in time, which means that some migmatites may crystallize at a different stage of anatexis than others (Urtson and Soesoo, 2009). This can lead to differences in how the migmatites look in field and considerably affect leucosome size measurements. Because of the abundance of melt-forming mineral triple junctions in many metapelitic rocks, a metapelite in an early stage of anatexis can contain relatively many small melt pockets (Sawyer, 2008). If the rock were to crystallize at this early stage of partial melting, the resulting leucosome width distribution would correspond with a truncated power-law distribution that has high exponent m and only spans a short size range ( Figure 6A).
Conversely, if the rock were to continue partially melting, the melt-bearing veins, and pockets would enlarge and merge . As depicted in Figure 6B, solidification at this stage would result in a leucosome size distribution with a lower exponent m and a wider size range than earlier, as there would be more mid-sized and large leucosomes than previously. This could appear either as a truncated or an untruncated power-law distribution, depending on if the maximum sizes of leucosomes in the system were of the observed (outcrop) scale or of a larger, unobservable scale. In continued partial melting, some of the largest melt-bearing veins or batches would grow large enough to escape the host rock . If the melt escape occurred for each melt-bearing vein in a certain maximum width range, e.g., decimeter scale, the resulting leucosome width distribution for a migmatite solidified after the removal of the melt would follow a truncated power-law. In this case, as Figure 6C illustrates, leucosomes smaller than the maximum size range would follow a power-law as expected, and the absence of large melt-bearing veins would be indicated by the truncation at the maximum size range widths.
If, on the other hand, there was a specific size range that melt accumulations generally grew up to before many of them escape at a single melt loss event, a migmatite solidified right before the melt loss event could occur would display a positive kink in leucosome size distribution. Figure 6D shows that the leucosomes around the kink-size would seem overrepresented in comparison to both smaller and larger leucosomes. In contrast to a truncated power-law trend, a power-law like distribution can be identified from leucosomes on either side of the kink, indicating that the two size ranges (melt accumulations smaller than the kink vs. those larger than the kink) developed without continuous interaction. According to this model, double power-law distributions with negative kinks could display the state of the system after kink-sized accumulations were emptied, as illustrated in Figure 6E.
It is unclear what could cause melt entrapment in particularly sized conduits. As earlier leucosome size distribution research in other areas has shown single power-law distributions rather than double power-law distributions, the double power-law distributions observed here may stem from some intrinsic properties of the Olkiluoto migmatite that are not present in all migmatites. Further research is certainly needed to ascertain how common double power-law distributions are in different migmatite areas and how the presence of kinks correlates with the physical properties of paleosomes and assumed protoliths. For example, the roles of paleosome grain size and structural anisotropy of the host rock in imposing physical constraints to melt mobility should be studied, as well as the effect of reversals relative strength of melts, leucosomes, and hosts during the whole anatectic process. Coupled with deformation, water diffusing away from melt can increase the efficiency of both crystallization and residual melt expulsion (White and Powell, 2010), and further research could establish if this process can lead to leucosome crystallization preferring certain sizes. As shear zones commonly occur in anatectic terranes, and they have also been shown to follow power-law distributions (de Riese et al., 2019), the interplay between the self-organization of strain and that of anatectic melt is another line of research to explore. Furthermore, identification of possible kinks and multiple power-law domains in other size ranges of the anatectic processes, for example in pluton scale as in the work by Soesoo and Bons (2015), would be necessary to understand which other processes can, at least temporarily, halt the self-organization of the anatectic system.

CONCLUSION
Because migmatites always represent a range of snapshots of the anatectic processes from the time of solidification, anatexis may be arrested at a stage where it does not appear to follow a single power-law distribution despite being a self-organized process. This can be due to temporal differences so that one self-organized anatectic system, with leucosome sizes following a single powerlaw distribution, solidified before another self-organized anatectic system was active in the same host rock and caused a sort of overprint of its leucosomes on top of the earlier system. Alternatively, double power-law distribution of leucosomes can depend on some disturbance in the self-organization of the anatectic system, causing melts to stall in certain-sized conduits before they can be emptied at once.

DATA AVAILABILITY STATEMENT
The leucosome width data measured and analyzed in this study can be found in Mendeley doi:10.17632/czchjrwrn2.1.

AUTHOR CONTRIBUTIONS
Conceptualization by ASa and ASo; methodology by ASa, ASo and CA; formal analysis by ASa and CA; investigation and visualization by ASa; writing of the original draft by ASa, CA and KN; review and editing by all authors; supervision by ASo, KN and OE; funding acquisition by ASa and OE. All authors contributed to the article and approved the submitted version.