Laurent Guillou-Frottier

Evgueni Burov

Thierry Augé

Eric Gloaguen

Rheological conditions for emplacement of Ural-Alaskan-type ultramafic complexes

Keywords: Ural-Alaskan complexes, Dunite core geometry, Magma ascent mechanisms, Platinum-group elements, Numerical modelling, Rheology, Softening

Ural-Alaskan-(or Alaskan-) type complexes correspond to a particular class of ultramafic intrusions that attract particular attention due to their deep mantle origin and their platinumgroup element (PGE) mineralization. When defined as massifs of dunite-clinopyroxenite, only forty-six complexes are reported in the literature. These large-scale dunite pipe-like structures are rarely isolated and they even can appear in clusters. To better understand genesis of these relatively young (< 460 Ma) complexes, a worldwide compilation has been built, and three categories have been defined: single circular or elliptical bodies, twin bodies with similar shapes, and dismembered dunite bodies. PGE enrichment in Alaskan-type complexes is highest for the second category, where twin bodies are interpreted as horizontal sections of Y-shaped dunite pipes. To constrain mechanical properties of the lithosphere allowing emplacement of the Alaskan-type complexes, the forceful diapiric ascent hypothesis is investigated through numerical thermo-mechanical models. One hundred high resolution experiments accounting for realistic phase changes and softening mechanisms have been performed. The experiments show that with no rheological softening of the host rock and in case of a relatively weak ductile lower crust, the uprising magma tends to spread laterally without reaching the surface. To account for the forceful ascent of deep magmas, it is hence necessary to assume a strong lower crust rheology and strong local softening mechanisms. Besides reproducing the clustered distribution of the weakness zones representing magma pathways, these latter experiments reproduce large-scale pipe-like (cylindrical) structures, Y-shaped and funnel-shaped bodies, and laterally-shifted structures.

Interestingly, zones of highest strain rates are located at the bottom parts of the inclined edges of Y-shaped and funnel-shaped bodies. The restricted age range of Alaskan-type complexes (< 460 Ma) would mean that prior to this time, the lower crust was less resistant due to the hotter geotherm, prohibiting the possibility of "Alaskan-type magmatism".

A C C E P T E D

M A N U S C R I P T

Introduction

Alaskan-type (or Ural-Alaskan-type) complexes represent a group of ultramafic bodies where orthopyroxene is generally absent, and can thus be described as duniteclinopyroxenite massifs. They have been originally defined as intrusions of olivine-rich lithologies (considered as mantle-derived material) with a concentrically zoned distribution, evolving from a dunite core to a clinopyroxenite rim, with gabbros or diorites at the margins [START_REF] Taylor | The zoned ultramafic complexes of southeastern Alaska, Part 4[END_REF]. The size of Alaskan-type complexes ranges from a few to hundreds of km², with one exception, the Guli massif, reaching 2000 km². In map view, they present circular or elliptical geometries, and can appear as isolated bodies or accompanied by a twin body of similar size, each of them being separated by a few kilometres. At a larger scale (hundreds of kilometres), several Alaskan-type bodies can be gathered in clusters, like in the Urals or in south-eastern Alaska (Figure 1). Although ultramafic rocks of various typologies can be found in many places on Earth, Alaskan-type complexes are very scarce since only 46 complexes can really fit the above definition of a dunite core surrounded by a clinopyroxenite rim (Table 1). Twelve other ultramafic bodies were classified in the past as Alaskan-type complexes but they do not include a dunite core, and hence they are not listed in Table 1 (e.g. Turn Mountain, Salt Chuck and Sukkwan Island ultramafic bodies in SE Alaska, [START_REF] Himmelberg | Characteristics and petrogenesis of Alaskan-type ultramaficmafic intrusions, southeastern Alaska[END_REF]. There is also no evidence for a dunite core in the so-called Alaskan-type complexes of Gnat Lake or Menard Creek ultramafic bodies in British Columbia [START_REF] Nixon | Geology and platinum-groupelement mineralization of Alaskan-type ultramafic complexes in British Columbia[END_REF]. Similarly, although [START_REF] Pettigrew | The Quetico intrusions of Western Superior Province: Neo-Archean examples of Alaskan/Ural-type mafic-ultramafic intrusions[END_REF] suggested that the Neo-Archean ultramafic intrusions in the Quetico subprovince, Canada, correspond to Alaskantype complexes, it appears that their cores are made of wehrlite (less than 90% of olivine)

and not of dunite. Nevertheless, even if the listed number of Alaskan-type complexes depends on their definition and on the interpretation of geological maps, their scarcity is noteworthily evident. Besides their scarcity, dunite-clinopyroxenite massifs show another intriguing feature, dealing with their ages: Alaskan-type complexes were emplaced during a small part of the Earth history, from 460 to 20 Ma.

A C C E P T E D M A N U S C R I P T

Ultramafic bodies and Alaskan-type complexes in particular, have been studied for a long time because of their associated mineralization in platinum-group elements (PGE) [START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF]. Platinum-Group Elements, which refer to six metals (osmium, iridium, ruthenium, rhodium, platinum and palladium), have become essential to modern industry and are now classified as "strategic" metals. Among these elements, two of them, iridium and platinum, are particularly present in Alaskan-type bodies, and it appears that the price of these two metals has been continuously increasing for the last decade. However, the number of exploited PGE ore deposits is largely limited to South Africa (Bushveld complex) and Russia (Noril'sk intrusion). In addition, the principal source of PGE in Alaskan-type complexes is placer deposit and not primary mineralization. This partly explains that there are only few studies dedicated to the emplacement mechanisms of Alaskan-type complexes which could help understand the conditions of PGE mineralization.

Platinum "nuggets" have been discovered and exploited within placers spatially associated with Alaskan-type bodies (e.g. [START_REF] Slansky | Platinum mineralization in the Alaskan-type intrusive complexes near Fifield, N.S.W., Australia. Part 2. Platinum-group minerals in placer deposits at Fifield[END_REF][START_REF] Tolstykh | The association of platinum-group minerals in placers of the Pustaya River, Kamchatka, Russia[END_REF]Malitch et al., 2002;[START_REF] Barkov | Platinum-group minerals from five placer deposits in British Columbia, Canada[END_REF]. Placers remain an important source for PGE production in Russia, after the Ni-Cu sulfides ores of Noril'sk-type intrusions [START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF][START_REF] Malitch | Pt-Fe nuggets derived from clinopyroxenite-dunite massifs, Russia: a structural, compositional and osmium-isotope study[END_REF]. According to [START_REF] Thakurta | Geochemical constraints on the origin of sulphide mineralization in the Duke Island complex, southeastern Alaska[END_REF] who studied the Alaskan-type Duke Island complex, southeastern Alaska, "the potential for the occurrence of world-class Cu-Ni-PGE sulphide deposits in Alaskan-type complexes should not be overlooked". Although some authors consider the absence of sulphides as representative of Alaskan-type bodies (e.g. [START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF], this feature is not retained in our definition (a dunite core surrounded by clinopyroxenites).

Another recent study [START_REF] Nazimova | The Galmoenan massif, Kamchatka, Russia: geology, PGE mineralization, applied mineralogy and beneficiation[END_REF] emphasized on the important platinum mineralization of the Galmoenan massif, Kamchatka, possibly defining, in addition to its wellknown placer deposits, an economic, medium-size platinum deposit of primary type. Two hundred kilometres southwest of the Galmoenan complex, numerous platinum "nuggets" were also found in placers from the Pustaya River, which crosses the Koryak-Kamchatka platinum belt related to eroded Alaskan-type intrusions [START_REF] Tolstykh | The association of platinum-group minerals in placers of the Pustaya River, Kamchatka, Russia[END_REF]. It must be

A C C E P T E D M A N U S C R I P T

stressed that placers around Alaskan-type bodies are not the sole type of PGE mineralization. Indeed, several examples have shown that primary type mineralization, where PGE are associated with chromitite lenses within dunites, have been exploited in the past (e.g. Nizhny Tagil, [START_REF] Augé | Primary platinum mineralization in the Nizhny Tagil and Kachkanar ultramafic complexes, Urals, Russia: a genetic model for PGE concentration in chromite-rich zones[END_REF], but their contribution remains historically very low compared to placer production. However, chromite-bearing dunites within restricted parts of the dunite cores appear to represent a potential economic target (e.g. Galmoenan, Nazimova et al., 2011).

Often associated with the location of ancient subduction zones (Urals, Alaska, British Columbia, Kamchatka, Colombia), Alaskan-type complexes can also form in intracontinental settings (Siberia, Ethiopia, Far East Russia). Several mechanisms have been invoked to explain the ascent and emplacement of mantle-derived magmas at the Earth's surface.

Melting from mantle plume heads or from subduction settings have both been proposed.

Although Alaskan-type complexes are not necessarily linked to subduction, they are still considered by some researchers as geological imprints of past subduction zones (e.g. [START_REF] Su | Occurrence of an Alaskan-type complex in the middle Tianshan massif, Central Asian Orogenic belt: inferences from petrological and mineralogical studies[END_REF]. On the contrary, Sutherland (1998) suggested a plume model where magmatic diapirs derived from the plume head rise in zones of increasing tensional stress, creating "tube-like" intrusions (see also [START_REF] Burg | Infra-arc mantle-crust transition and intra-arc mantle diapirs in the Kohistan Complex (Pakistani Himalaya): petrostructural evidence[END_REF][START_REF] Pirajno | The association of mafic-ultramafic intrusions and A-type magmatism in the Tian-Shan and Altay orogens, NW China: implications for geodynamic evolution and potential for the discovery of new ore deposits[END_REF]. Studying ascent mechanisms of ultramafic bodies, [START_REF] Gerya | Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation[END_REF] and [START_REF] Burg | Translithospheric Mantle Diapirism: Geological Evidence and Numerical Modelling of the Kondyor Zoned Ultramafic Complex (Russian Far-East)[END_REF] found that large-scale geodynamic processes are not necessarily required, and that local rheological properties could lead to "translithospheric diapirism", a mechanism with which mantle-derived material could quickly rise up through the lithospheric mantle and the crust. In particular, a cold elasto-plastic crust would facilitate the upward propagation through the crust. This process was initially suggested by [START_REF] Bodinier | The Konder « Anorogenic » Ultramafic Massif (Aldan Shield) : Impact of a Translithospheric Mantle Diapir at the Earth's Surface ? Fourth International Workshop on Orogenic Lherzolite and Mantle Processes[END_REF], who analysed microstructural, geochemical and mineralogical data from dunite and pyroxenites of the Kondyor complex.

The authors concluded from previous gravimetric data and tectonic reconstructions that the massif represents "an apex of a dense ultramafic cylinder rooted down to, at least, 10 km depth". Magma at the origin of the Kondyor dunite would have accumulated at the base of the lithospheric mantle, triggering a forcefully ascent of a solid-state dunite diapir. Here and

A C C E P T E D M A N U S C R I P T
in the following, the "solid-state" ascent refers to the rheological behaviour of mantle-derived material, which corresponds to partially molten peridotite with a low density but a high effective viscosity.

The present study consists in looking for appropriate rheological conditions allowing simulations of magmatic processes at the origin of Alaskan-type complexes (fast ascent through the entire lithosphere, particular geometries at the surface, clustering process), through lithospheric-scale thermo-mechanical models that account for the properties of the involved magma. Realistic rheological laws for the crust and the lithosphere (see Supplementary material), as well as realistic phase changes, have been included in the thermo-mechanical numerical model. In particular, our numerical results focus on spatial distribution and geometries of the weakness zones that are created by the ascent of a deep and hot mafic magma.

Alaskan-type complexes: general description and classification

Alaskan-type complexes are defined by the following successive events: : first a dunite core is emplaced; second, melt-rock reactions develop successively rims of wehrlite, followed by clinopyroxenites ± hornblendite and then monzonite-gabbro (i.e. absence of orthopyroxene); the contacts between ultramafic and surrounding gabbroic rocks are sharp;

the chemical composition of olivine in dunites is highly magnesian; chromite is a phase occurring exclusively in dunite, where it is ubiquitous and locally forms concentration occurring as schlieren or pseudo-layers enriched in PGE (see [START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF] for a detailed review). Ultramafic and mafic rocks contain significant proportions of hydrated minerals (amphibole, phlogopite) and hydrated bodies (ultramafic-mafic pegmatoids). Additional common features to Alaskan-type complexes were later suggested, like the specific concentration and distribution of rare-earth elements (REE) in ultramafic and gabbroic rocks, or the relative absence of sulphides (see [START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF]. In the following, we present a distinct classification, based on pure geometrical considerations.

A C C E P T E D M A N U S C R I P T

When geological maps of Alaskan-type complexes are examined, three distinct classes based on surface geometries can be defined, regardless of the geodynamic setting (Figure 2). The first category includes circular or elliptical surface expression of dunite and clinopyroxenite rims with an external diameter of several kilometres (Figure 2a). For some cases, more elongated shapes are observed (e.g. Yubdo, Gaositai, Nizhny Tagil, Wrede Creek, Figure 2c). The second category corresponds to twin bodies (similar shape and size),

showing dunite cores of several kilometres in size separated by a few kilometres, each of them being surrounded by clinopyroxenites (e.g. Galmoenan, Goodnews Bay and Duke Island in Figure 2b; Kachkanar, Owendale and Condoto in Figure 2d). The third category includes disconnected dunite blocks of different sizes (e.g. Lunar Creek and Uktus, Figure 2d). All the other Alaskan-type complexes not shown in Figure 2 but listed in Table 1 can be classified in one of the three categories. For example, geological maps of the Turnagain, Red

Bluff or Kane Peak ultramafic bodies [START_REF] Himmelberg | Characteristics and petrogenesis of Alaskan-type ultramaficmafic intrusions, southeastern Alaska[END_REF][START_REF] Scheel | Chromian spinel in the Turnagain Alaskan-type ultramafic intrusion, northern British Columbia, Canada[END_REF] clearly indicate that they belong to the first category of Alaskan-type complexes. In British Columbia, the Polaris complex belongs to the second category since it comprises two main dunite cores separated by wehrlite and clinopyroxenite [START_REF] Nixon | Geology and platinum-groupelement mineralization of Alaskan-type ultramafic complexes in British Columbia[END_REF]. Third column of Table 1 indicates which category the Alaskan-type complex belongs to, with the number of dunite cores (one, two or several), when available.

Vertical models

As mentioned above, gravity measurements around the Kondyor massif led Efimov and Tavlin (1978) to suggest that ultramafic rocks extend as a cylindrical pipe down to a depth of at least 10 km. Interpretations of geological data at Union Bay, southeastern Alaska, by [START_REF] Ruckmick | Origin of the ultramafic complex at Union Bay, south-eastern Alaska[END_REF] and [START_REF] Himmelberg | Characteristics and petrogenesis of Alaskan-type ultramaficmafic intrusions, southeastern Alaska[END_REF] suggest that ultramafic rocks delineate an almost vertical pipe with slightly decreasing diameter with depth (funnelshaped pipe). Again, based on geological observations, van der Poel and Hinderman (2000)

described the two twin bodies of the Goodnews Bay complex as vertical pipes widening out at the surface. Recently, the geometric model of the Duke Island complex proposed by [START_REF] Thakurta | Geochemical constraints on the origin of sulphide mineralization in the Duke Island complex, southeastern Alaska[END_REF], which was derived from a previous model by [START_REF] Butler | Paleomagnetism of the Duke Island, Alaska, ultramafic complex revisited[END_REF],

A C C E P T E D M A N U S C R I P T
illustrates a three-branches structure, two of which including dunite pipes surrounded by clinopyroxenites, joining at depth and forming a common "trunk". To summarize, vertical structure of Alaskan-type complexes would correspond to magmatic conduits (pipes) which could have pierced the mantle lithosphere and the crust through a diapirism mechanism [START_REF] Bodinier | The Konder « Anorogenic » Ultramafic Massif (Aldan Shield) : Impact of a Translithospheric Mantle Diapir at the Earth's Surface ? Fourth International Workshop on Orogenic Lherzolite and Mantle Processes[END_REF][START_REF] Gerya | Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation[END_REF][START_REF] Burg | Translithospheric Mantle Diapirism: Geological Evidence and Numerical Modelling of the Kondyor Zoned Ultramafic Complex (Russian Far-East)[END_REF].

Clusters

Clustering of Alaskan-type complexes (Figure 1b) can be quantified by the averaged distance between nearest neighbours, which is around 60 km for the Urals, and around 30 km for south-eastern Alaska. For both cases, the role of past or present subduction zones has been invoked [START_REF] Murray | Zoned ultramafic complexes of the Alaskan type: feeder pipes of andesitic volcanoes[END_REF]. Actually, there is no reason to get same values since volcanism related to subduction seems to be randomly distributed [START_REF] Brémond D'ars | Distribution of volcanoes in active margins[END_REF]. According to these authors, spatial distribution of subduction-related volcanoes would be representative of melting conditions and deep magmatic activity, such as recurrent magmatic pulses.

Clustering of Alaskan-type complexes is not necessarily linear (like in the Urals, Figure 1) and several neighbouring bodies seem to have been emplaced in a single event (see in Table 1 the ages of Alaskan-type complexes in the Urals, in Kamchatka, and in SE Alaska). Simultaneous emplacement of such deep magmatic bodies probably requires a particular thermo-mechanical event affecting the bulk crust over several hundreds of km.

Uprising of mantle-derived material through the lithosphere requires the development of weak crustal zones. The clustering phenomenon (of weak zones) might be analogous to the faulting spatial periodicity problem. Fracturing and faulting spatial periodicity has been extensively studied in the past decades (e.g. [START_REF] Ouillon | Hierarchical geometry of faulting[END_REF]. These authors concluded that faulting spatial periodicity is mainly controlled by rheological layering of the crust and not by the complex rupture process. Similarly, it is the rheology of the lower crust which appears to control folding of the lithosphere [START_REF] Gerbault | Prediction of faulting from the theories of elasticity and plasticity: what are the limits ?[END_REF]. In the following sections, the clustering of Alaskan-type complexes is investigated through distinct possible crustal rheologies.

A C C E P T E D M A N U S C R I P T

Physical properties of ultramafic rocks and magmas

Physical properties of olivine-rich ultramafic rocks (density, thermal properties, and rheology) have been studied by several teams in the past [START_REF] Chopra | High-temperature transient creep in olivine rocks[END_REF]Paterson, 1981, 1984;[START_REF] Chopra | High-temperature transient creep in olivine rocks[END_REF]Mei andKohlstedt, 2000a, 2000b;[START_REF] Suzuki | Density of peridodite melts at high pressure[END_REF][START_REF] Gibert | Thermal diffusivity of olivine single crystals and a dunite at high temperature: evidence for heat transfer by radiation in the upper mantle[END_REF][START_REF] Wang | A review of water contents and ductile deformation mechanisms of olivine: implications for the lithosphere-asthenosphere boundary of continents[END_REF][START_REF] Keefner | Dependence of dislocation creep of dunite on oxygen fugacity: implications for viscosity variations in Earth's mantle[END_REF]. Recently obtained data on clinopyroxene crystals from the Konder and Inagli massifs have shown a close composition and a similar configuration of REE patterns to those obtained with samples from the Alaskan-type Kytlym and Nizhni Tagil massifs [START_REF] Simonov | Genesis of platiniferous massifs in the southeastern Siberian platform[END_REF]. The authors concluded that dunites of the Konder and Inagli massifs were produced by water-rich (up to 0.6%) magmas at temperatures of 1300-1460°C. This is in accordance with the amount of hydrated minerals (amphibole, phlogopite) and hydrated bodies (ultramafic-mafic pegmatoids) observed in such Alaskan-type complexes.

It is well-known that the presence of water strongly modifies rheological properties, and such hydrated magmas at high temperature may indeed undergo unexpected flow regimes. The presence of water in olivine rich magmas also induces a decrease in viscosity:

a reduction factor of 4 to 10 should be applied to mantle-derived material if viscosity of midoceanic ridge basalts had to be used (Mei and Kohlstedt, 2000b). Laboratory experiments attempting to reproduce realistic thermal and mechanical conditions are however limited. As emphasized by [START_REF] Korenaga | A new analysis of experimental data on olivine rheology[END_REF], olivine rheology is complex since a minimum of eighteen parameters are required to write rheological laws. In addition, rheological laws for water-rich olivine suffer major uncertainties, as it is actually the case for crustal rocks [START_REF] Burov | Plate rheology and mechanics[END_REF]. Recently, high-pressure and high-temperature experiments on deformation of dunite demonstrated that accounting for the role of oxygen fugacity in the range of upper mantle conditions, could decrease dunite viscosity by a factor of ~5 [START_REF] Keefner | Dependence of dislocation creep of dunite on oxygen fugacity: implications for viscosity variations in Earth's mantle[END_REF]. [START_REF] Korenaga | A new analysis of experimental data on olivine rheology[END_REF] illustrate different effective viscosity profiles for the upper 300 km of the mantle, for dry and wet cases and for different chosen geotherms; all viscosity profiles for

A C C E P T E D M A N U S C R I P T
wet conditions show uncertainties of several orders of magnitude. The authors particularly insist on the large scatter of experimental data under wet conditions.

As far as thermal properties are concerned, it seems that they can be underestimated when experimental data are considered. [START_REF] Gibert | Thermal diffusivity of olivine single crystals and a dunite at high temperature: evidence for heat transfer by radiation in the upper mantle[END_REF] concluded from their measurements that thermal diffusivity of mantle-derived material could be 50% higher than previously assumed. It follows that heat transfer processes related to low viscous olivine-rich rocks at high temperatures could thus be much more efficient, and thus much faster than commonly thought.

All the above-mentioned experimental data suggest that olivine-rich mantle rocks at high temperature may exhibit particular rheological behaviour involving strong strain and/or strain-rate softening mechanisms (see also [START_REF] Burov | Rheology and strength of the lithosphere[END_REF]. As far as modelling parameters are concerned, density is here determined from mineralogical phases defined by ambient thermo-mechanical conditions (pressure and temperature) through equations of state derived for corresponding mineralogical compositions (see below). The effective viscosity is assessed from rock mechanics data on flow law parameters for given thermo-mechanical conditions, together with testing parameters of softening rules defining local strain-induced reduction in viscosity (see next section). Thus, the rheological model comprises distinct softening rules that are tested until a reasonable translithospheric diapirism mechanism is achieved.

Numerical model and preliminary tests

In the following, we present a series of numerical experiments aimed at reproducing the apparent "lithospheric diapirism" which is assumed to explain the ascent and emplacement mechanisms of Alaskan-type complexes. Main model features and tests experiments are first presented (Figure 3). Then, the results obtained with specific softening rules allowing the reproduction of typical geometries of Alaskan-type bodies are discussed.

Model features

A C C E P T E D M A N U S C R I P T

We basically follow the philosophy of the earlier study by [START_REF] Gerya | Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation[END_REF], yet for a wider range of thermo-tectonic conditions, rheological assumptions and configurations specific to Alaskan-type complexes. We also test different softening rules and other parameters in order to assess the rheological properties of the lithosphere enabling for magma ascent from great depth. On the technical side, the main difference of our approach from that of the one developed in [START_REF] Gerya | Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation[END_REF] refers to a formally higher degree of physical consistency achieved through direct implementation of free surface, viscouselastic-plastic rheology and phase changes (the two latter features were also integrated in the latest versions of the code by T. Gerya but were not applied to the problem of intrusion of ultramafic magmatic bodies). In particular, incorporation of continuous phase changes [START_REF] Burov | Continental plate collision, P-T-t-z conditions and unstable vs. stable plate dynamics : Insights from thermo-mechanical modelling[END_REF] implies that rock density is thermo-dynamically consistent being dynamically updated as function of P-T (pressure and temperature) conditions from the equations of state  = f(P,T) based on minimization of Gibbs free energy for given mineralogical composition. With that goal, the thermodynamic algorithm PERPLEX and its data base [START_REF] Connolly | Computation of phase equilibria by linear programming: a tool for geodynamic modelling and its application to subduction zone decarbonation[END_REF] have been coupled with the main code. This approach not only provides more reliable values for density and for hence thermodynamically consistent definitions of thermal expansion  (P,T) and compressibility  (P,T), but also for thermal and elastic parameters, the amount of released or absorbed fluids and other physical properties that can be used for thermo-mechanical calculations and on post-processing stage. As mentioned above, PERPLEX minimizes free Gibbs energy G for a given chemical composition to calculate an equilibrium mineralogical assemblage for given P-T conditions:

1 n ii i GN     (1)
where  i is the chemical potential and N i the moles number for each component i constitutive of the assemblage. Given the mineralogical composition, the computation of density is then straightforward (e.g., [START_REF] Yamato | Burial and exhumation in a subduction wedge: Mutual constraints from thermomechanical modeling and natural P-T-t data (Schistes Lustrés, western Alps)[END_REF][START_REF] Burov | Continental plate collision, P-T-t-z conditions and unstable vs. stable plate dynamics : Insights from thermo-mechanical modelling[END_REF]. Thermal expansion  eff is corrected for latent heat following [START_REF] Gerya | Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation[END_REF]:

A C C E P T E D M A N U S C R I P T 0 ( , ) L eff T QM PT Tp          (2)
where  0 , Q L , M, T and p are density, latent heat of melt, volume fraction of melt, temperature and pressure, respectively. To solve the thermo-mechanical part of the problem, we use the code Flamar 12 (outgrowth of FLAC and Paravoz [START_REF] Cundall | Numerical experiments on localization in frictional materials[END_REF][START_REF] Poliakov | An explicit inertial method for the simulation of visco-elastic flow: an evaluation of elastic effects on diapiric flow in two-or threelayers models[END_REF]).

Flamar is a well-tested code that have been used for many geodynamic applications including exhumation of metamorphic domes [START_REF] Tirel | Dynamics and structural development of metamorphic core complexes[END_REF], mantle-lithosphere interactions [START_REF] D'acremont | Numerical modelling of a mantle plume: the plume headlithosphere interaction in the formation of an oceanic large igneous province[END_REF][START_REF] Burov | Controls of mantle plumes and lithospheric folding on modes of intraplate continental tectonics: differences and similarities[END_REF], continental collision (e.g., [START_REF] Toussaint | Tectonic evolution of a continental collision zone: a thermo mechanical numerical model[END_REF][START_REF] Burov | Continental plate collision, P-T-t-z conditions and unstable vs. stable plate dynamics : Insights from thermo-mechanical modelling[END_REF] and metamorphic rock exhumation (e.g., [START_REF] Burov | Mechanisms of continental subduction and exhumation of HP and UHP rocks[END_REF]. As its predecessors, Flamar is a large strain fully explicit time-marching algorithm that locally solves full Newtonian equations of motion in a continuum mechanics approximation:

0      u σg (3)
coupled with viscous-elastic-plastic constitutive equations (see Supplementary material):

( , , , ,... ...)

D FT Dt    u u u (4)
with equations of heat transfer and with state equations:

 C p DT/Dt -(k T) -H = 0 (5)  = f (P,T) (6)
Here u, , g, k are the respective terms for displacement, stress, acceleration due to body forces and thermal conductivity; dot represents the time-derivative; pressure P is negative for compression. In the eq. 5, heat advection term

T  u
is included in the Lagrangian derivative DT/Dt. The triangular brackets in (3) specify conditional use of the related term as in quasistatic mode inertial terms are damped using inertial mass scaling [START_REF] Cundall | Numerical experiments on localization in frictional materials[END_REF]. The terms t, C p, H designate respectively time, specific heat, and internal heat production (radiogenic heating, shear heating, latent heat). The expression  = f(P,T) refers to the density computed

A C C E P T E D M A N U S C R I P T
by a thermodynamic module PERPLEX as function of P and T. The terms D/Dt and F are the objective Jaumann stress time derivative and a functional, respectively. In the Lagrangian method, incremental displacements are added to the grid coordinates allowing the mesh to move and deform with the material. This allows for the solution of large-strain problems while using locally the small-strain formulation: on each time step the solution is obtained in local coordinates, which are then updated in a large strain mode.

Solution of ( 3) provides velocities at mesh points used for computation of element strains and of heat advection

T  u
. These strains are used in (4) to calculate element stresses and equivalent forces used to compute velocities for the next time step. Due to the explicit approach, there are no convergence issues, which is rather common for implicit methods in case of non-linear rheologies. The algorithm automatically checks and adopts the internal time step using 0.1 -0.5 of Courant's criterion of stability, which warrants stable solution.

Model setup and test experiments

We first run a series of experiments roughly following the strategy chosen by Gerya and Burg ( 2007), yet for conditions and parameters compatible with the Alaskan-type complexes. The setup of the models is following:

In the reference model, the computational domain is 1000 km wide and 300 km deep with a homogeneous grid comprising 600 x 260 rectangular elements. Some experiments were run using local mesh refinement allowing for higher spatial resolution in center of the model (0.8 km x 0.5 km). Depending on the experiment, the lithosphere was either 100 km or 200 km thick including 40 km thick crust. Semi-spherical 200 km (in some experiments 50 km or 100 km) wide magmatic source is placed at the bottom of the model. An initially narrow -2 element widedike-like conduit connects the source region with the bottom of the crust (dashed white contour in Figure 3a). Density and viscosity contrasts impose large and fast modifications in the lithospheric and crustal stress regimes; thermo-mechanical evolution of the system is thus computed until high strain rates are achieved at the surface. Composition,

A C C E P T E D M A N U S C R I P T
stress, strain and strain rates are computed for each grid cell, but "magma flow" (in terms of fluid dynamics) is not simulated in these experiments due to their time scales. The conduit is initially blocked at Moho depth to investigate further initialization of the potential magma pathways in the overlying crust. We do not continue experiments once the pathways are initialized since further propagation of magma and its interaction with the walls of the pathways are beyond the scope of this study, requiring a separate parametric investigation at local scale. The following illustrations show compositional fields and strain rate values that are computed after thousands years of evolution. In particular, large values of strain rates are assumed to be located where contrasts in physical properties are the highests. In other words, the terms "magma ascent" in the following text correspond to the upward migration of high strain rate zones.

In the reference experiments (Figure 3b

The influence of compositional density contrast.

The density of partially molten magma is subject to large variations that depend on composition, transient P-T conditions and fluid content. We thus first tested end-member cases with various density contrast cut-offs (as defined by ( m - b )/ b , where  m is magma density and  b the background rock density): from 0.1 to 0.5. In these experiments we did not activate thermodynamic phase changes. These experiments show the extreme importance of the density contrast on magma emplacement and ascent rate (see Figure 3b, where density contrast between the magma and the mantle equals 0.4 and 0.2). It is important to note that magma ascent velocity is extremely high, reaching several tens to hundreds of m/yr. These values are in accordance with estimates from previous analytical and numerical studies that account for realistic (non-Newtonian) rheologies (e.g. Weingerg and Podladchikov, 1994;[START_REF] Van Keken | Evolution of starting mantle plumes: a comparison betwen numerical and laboratory models[END_REF][START_REF] Burov | Ascent and emplacement of magma reservoirs in brittle-ductile upper crust[END_REF].

Experiments with thermodynamically consistent densities

The next series of reference experiments have been implemented with thermodynamically consistent account for phase changes (Figure 3c). Due to P-T dependent density, surface impact of magma is moderate compared to the previous models with fixed density contrast; a relatively large secondary magma region forms at the base of the crust, sourcing magma propagation to the surface. One can observe, similarly to the experiments of [START_REF] Gerya | Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation[END_REF], that magma ascent creates a deviatoric stress field provoking brittle faulting at surface. The created faults go deep down to 40 km depth and thus may serve for channeling magma to the surface. It is noteworthy that formation of surface faults is conditioned by 2 major factors: (1) the assumed crustal geotherm and (2) magma viscosity cut-off. The background geotherm should be relatively cold in the crustal region to make brittle deformation possible while magma should not be neither too viscous nor too fluid: its viscosity must be lower than 10 18 Pa.s, otherwise its ascent will be too slow and magma would be stalled at some depth, but greater than 10 15 Pa.s, otherwise magma might spread laterally somewhere in the crust instead of pushing the brittle crust up. These high viscosity values are one order of magnitude higher than those used by [START_REF] Gerya | Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation[END_REF], where the associated melt fraction was around 10 to 30% vol. However, in case of the presence of external extensional field (tectonic extension), faults may be formed even in cases where the above conditions are not satisfied (see an example in [START_REF] Burov | Thermomechanical behavior of large ash flow calderas[END_REF].

A C C E P T E D M A N U S C R I P T

Experiments with softening mechanisms

As can be seen from the experiments of the previous sections, the reference thermorheological structure results in relatively weak ductile lower crust, so that magmatic material ascending from the chamber spreads laterally at some level within the ductile lower crustal layer resulting in distributed deformation at surface since it is denser than the normal crust. It follows that the magmatic conduit is destroyed in its crustal part. Without detailing all experiments of this sensitivity study, we have found that preservation of the magmatic channel requires strong lower crustal rheology ("diabase/granulite") and strong "aggressive" softening of the brittle and ductile material. We have also extended the brittle failure criterion by adding to the common Mohr-Coulomb criterion that operates under overall compression (i.e. negative pressure, P < 0), the parabolic Griffith's criterion applicable in pure tension (positive pressure, P > 0):

 2 = 4T 0 2 +4T 0  n ( 7
)
where  is shear stress, T 0 is tension cut-off, typically 20-40 MPa,  n is normal stress. The aggressive softening for the brittle and brittle rock assuming similar relations to show that the observed morphologies of rifted basins cannot be explained without frictional-plastic and viscous strain softening. Similarly, [START_REF] Gerya | Intrusion of ultramafic magmatic bodies into the continental crust: numerical simulation[END_REF] and most other workers assumed softening laws that basically do not stem from experimental rock mechanics but from field observations and thermo-mechanical models that show that at macroscopic scale faults and shear bands should undergo strong softening to explain the observed strain patterns and evolution of geological structures. In general, linear strain-softening laws (e.g., [START_REF] Huismans | Asymmetric lithospheric extension: Relative importance of frictional-plastic and viscous strain softening inferred from numerical experiments[END_REF] are used and we here follow the same approach by implementing a simple piecewise linear ductile softening rule:

if    0 , () = (T,,  ' … ) (8) if  1 >  >  0 , () =  (1 -(- 0 )/(  1 - 0 ) ) (9) if    1 , () =  (1-) (10)
where  0 and  1 are respectively minimal and maximal accumulated strain values, within which progressive softening is applied;  is effective initial viscosity derived as function of T, ,  ' and material parameters;  is a constant viscosity factor such that (1-) is maximal ratio between viscosity after and before softening, e.g.,  = 0.99900 means 3 order maximal

A C C E P T E D M A N U S C R I P T

reduction of viscosity. We changed  0 from 0.1 to 0.2,  1 from 0.5 to 2.0 and  from 0.90000 to 0.99999. Seven particularly remarkable cases are detailed below.

Case 1.

The experiment shown in Figure 4 corresponds to  0 = 0.2,  1 = 0.5,  = 0.99900. This experiment clearly shows that the ascending magma has widened the initially thin channel (Figure 4a). Before arriving at the base of the crust, fracturesor weakness zonesdeveloped at several locations, as it can be seen in the strain rate pattern of Figure 4b, which also shows localization of ductile strain near channel walls that are mechanically eroded during magma ascent. Note in Figure 4c that the surface uplift extends over a distance greater than 100 km, and that spacing between crustal channels of high strain rates equals, respectively, 42, 23, 10, 20, 15 and 15 km (av. 20 ± 10 km).

An additional peculiar feature is revealed in Figure 4b: some zones of high strain rate values in the crust (light green and yellow zones) exhibit a triangular configuration, or funnelshape geometry. In addition, over the 8 fractured (weakness) zones, only one appears with an angle (laterally shifted body at shallow depths), the others being mostly vertical.

Interestingly, 4 of them are not directly connected to the magmatic channel. These observations are consistent with several characteristics of the previously described Alaskantype complexes, like the clustering phenomenon (see also discussion). Figure 4c illustrates the topography evolution, where fractures/channels are characterized by local subsidence, which is not visible at early stages. The large-scale doming begins to appear at ~0.5 kyr and clearly develops at later stages (after 2 kyr) while ascending magma has settled in the vertical conduits/channels and has sufficiently cooled.

Case 2.

These experiments (Figure 5) correspond to another softening mechanism, where controlling parameters are chosen as follows:  0 = 0.1,  1 = 1.0,  = 0.99900. Reduction in viscosity is thus similar to the previous experiment, but the range of the accumulated strain values, for which progressive softening is applied, becomes larger. In other words, the upper

A C C E P T E D M A N U S C R I P T

crust becomes weaker with time than in the previous experiment. It follows that the surface uplift observed before is smaller in these experiments, as can be seen from Figure 5.

However, pattern of fracturation, which is actually determined by dynamics of ascending magma, remains relatively similar, consisting in regularly spaced channels of high strain rates (Figure 5). In that case, the funnel-shape structure is not observed, but two Y-shaped structures are present at both edges of the channel tip. These Y-shaped channels differ from conjugate fault zones in the crust, since they could be characterized by larger distance between two branches, by larger vertical size and potentially smaller width, and could also be called, in that case of Figure 5, "U-shaped" structures. It can also be mentioned that topography shows a slight subsidence at the top of each channel, resulting in an undulated surface topography over 100-150 km.

Case 3 and 4.

The experiments shown in Figures 6a and6b correspond to another softening mechanism, where controlling parameters are chosen as follows:  0 = 0.2,  1 = 2.0,  = 0. 99999 (meaning that the reduction in viscosity reaches a factor of 10 5 ). The range of strains where softening is applied is thus ~ twice wider than that in previous experiments. This results in a formation of a volcanic-like structure, where deep magma erupts at the surface.

Figure 6a corresponds to a tension cut-off value of 10 MPa and Figure 6b shows the case of a tension cut-off value of 1000 MPa. Fractured zones can be inferred from step-like aspects of crustal layers.

Case 5.

In Figure 6c, controlling parameters are  0 = 0.2,  1 = 0.75,  = 0. 99900 and tension cut-off equals 1000 MPa (no Griffith's criterion). Here again, fracture zones can be easily seen with step-like structures in crustal layers (see inset of Figure 6c). Magma reaches the upper crust and fracturing extends over ~100 km.

The three cases shown in Figure 6 do not illustrate Alaskan-type structures since at the lower crust level, the width of the magmatic channel has been enlarged up to 30-40 km.

Contrarily to case 2 (shown in Figure 4) there is, for these 3 cases, no mechanical barrier

A C C E P T E D M A N U S C R I P T

above which high strain rate zones would become distributed at the surface over 100-150 km. The following tests show additionnal Alaskan-type structures which were obtained while controlling parameters ( 0 ,  1 ,  and T 0 ) were varied.

Case 6 and 7

Figure 7 shows two additional experiments where controlling parameters are respectively:  0 = 0.2,  1 = 1.0,  = 0.99900, and  0 = 0.1,  1 = 1.0 ,  = 0.99990. The first experiment (Figure 7a) show seven Y-shapes structures, whose "trunks" are separated by 28, 55, 25, 23, 30 and 55 km (av. 36 ± 13 km). The second case (Figure 7b) shows seven cylindrical to funnel-shaped structures, separated by 25, 28, 10, 20, 30 and 32 km (av. 24 ± 7 km). As for cases 1 and 2 (Figures 4 and5), a single deep magmatic source creates several weakness zones in the upper crust, separated by ~20-40 km, thus simulating the clustering aspect of magmatic channels (linear clusters shown in Figure 1).

Contrarily to cases 3, 4 and 5 (Figure 6), where deep material reaches the surface with no apparent mechanical barrier, the other results (Figures 4, 5 and7) clearly shows at ~20 km depth some drastic changes in the mechanical regime, that favour lateral distribution of high strain rate values.

Individual bodies

In addition to the reproduction of clustered structuresat least in these 2D numerical modelsthe obtained geometries of individual weakness zones (or areas with high strain rate values) appear to replicate the schematic models inferred from geological studies, such as the laterally-shifted bodies suggested for Union Bay and Goodnews Bay complexes [START_REF] Himmelberg | Characteristics and petrogenesis of Alaskan-type ultramaficmafic intrusions, southeastern Alaska[END_REF][START_REF] Van Der Poel | Summary report on lode platinum exploration[END_REF], or the Y-shaped body suggested by [START_REF] Thakurta | Geochemical constraints on the origin of sulphide mineralization in the Duke Island complex, southeastern Alaska[END_REF] for the Duke Island complex. Figure 8 illustrates distinct numerical results of such individual geometries, where it can be seen that dimensions (largely exceeding the resolution limit) correspond to natural cases (see Table 1). Moreover, the amplitudes of the strain rate values (from light green to red delineated by black dashed lines) provide additional information: (1) the strain rate values tend to be higher in funnelshaped and shallow Y-shaped structures; (2) these zones are located near the bottom parts

A C C E P T E D M A N U S C R I P T
of inclined edges. When there are no inclined edges at shallow (~5-10 km) depths, i.e. for vertical pipes such as the two left ones of Y-(or more exactly U-) shaped structures (middle cases in Figure 8), strain rate values are rather homogeneous within the entire structure.

Alaskan-type complexes and their PGE mineralization

6.1 PGE concentration PGE mineralization in Alaskan-type complexes is mainly associated with chromitites hosted by the dunite core. When compared to dunites, chromite concentrations appear to be significantly enriched in PGE by a factor of ~5-6 [START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF]. As emphasized by many studies, mineralized chromitites from Alaskan-type complexes exhibit PGE concentrations anomalously enriched in iridium and platinum (Figure 9). Indeed, the available data (with few exceptions discussed hereafter) show a typical "M" pattern, the two peaks representing Ir and Pt concentrations [START_REF] Garuti | Platinum-group elements as petrological indicators in mafic-ultramafic complexes of the central and southern Urals: preliminary results[END_REF][START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF][START_REF] Augé | Primary platinum mineralization in the Nizhny Tagil and Kachkanar ultramafic complexes, Urals, Russia: a genetic model for PGE concentration in chromite-rich zones[END_REF]Zacarini et al., 2011;[START_REF] Nazimova | The Galmoenan massif, Kamchatka, Russia: geology, PGE mineralization, applied mineralogy and beneficiation[END_REF], whereas, by comparison, mineralized podiform chromitites from ophiolites show a decreasing trend from Os-Ir-Ru to Pt-Pd elements (e.g. [START_REF] Ahmed | Unexpectedly high-PGE chromitites from the deeper mantle section of the northern Oman ophiolite and its tectonic implications[END_REF][START_REF] Büchl | Formation of podiform chromitite deposits: implications from PGE abundances and Os isotopic compositions of chromites from the Troodos complex[END_REF][START_REF] Zhou | REE and PGE geochemical constraints on the formation of dunites in the Luobusa ophiolite, southern Tibet[END_REF].

The PGE spidergram of Figure 9a illustrates this "M" pattern, as defined by typical mineralized chromites from Alaskan-type bodies. Spidergram of Figure 9b gathers other data which can be considered as coming from anomalously mineralized chromites (either an excess in Os or Pd, or an anomalously low signature in Pt). The Nizhny Tagil 3 and 4 samples correspond to chromitites of barren type and brecciated chromitites, respectively [START_REF] Augé | Primary platinum mineralization in the Nizhny Tagil and Kachkanar ultramafic complexes, Urals, Russia: a genetic model for PGE concentration in chromite-rich zones[END_REF]. Sample from the Butyrin vein in the Kytlym complex shows high concentrations of Os and Pd, while the two other samples (Kytlym 2 and 3) exhibit only slight anomalous Os concentrations and are thus not considered as anomalously mineralized. It is interesting to note that the distribution from the largest Alaskan-type body (the Guli massif, whose total area is around 2000 km², including two dunite cores covering 450 km²), exhibit the lowest Pt concentration with relatively high values for Os, Ir and Ru (Figure 9b). Malitch et al. (2002) also underlined that the common occurrence of Os-Ir-Ru alloys in chromitites

A C C E P T E D M A N U S C R I P T
from the Guli massif matches with patterns obtained with ophiolitic massifs. This fact suggests that body size and mineralization features may be related. Another similar observation is described below.

Using geometrical classification described in section 2, the above PGE concentrations from typical mineralized chromites of Alaskan-type bodies (Figure 9a) have been separated in two groups defined by the number of dunite cores (Figure 9c and9d). It turns out that single dunite cores (Figure 9c) show the lowest ranges of PGE concentrations (curves are located in the lower part of the M pattern area). Figure 9d gathers data from only twin Alaskan-type bodies, and include all the highest PGE concentrations. Note that data from the Uktus 2 sample, which may belong to Figures 9c or 9d (due to uncertainty in the initial number of dunite cores), do not alter this observation.

Inclined edges

Bottom parts of funnel-shaped ultramafic bodies (not necessarily of Alaskan-type) often correspond to mineralized areas (see discussion). Indeed, these high strain rate zones may concentrate magma pathways, thus enhancing potential for mineral deposition. These high strain rate zones (see section 5.6 and Figure 8), more frequent at the outer edges of the structure, may also explain the concentrically character of Alaskan-type complexes: in case of multiple magmatic pulses, the ascent of dunite cores would be easily followed by subsequent emplacement of pyroxenites around the dunite core, where high strain rate localized during the first ascent. Some field studies have shown a syn-emplacement deformation localized within the complex periphery (e.g. pyroxenite mylonites associated with dykes along borders of Svetloborsky complex, Urals, [START_REF] Tolstykh | Platinum mineralization of the Svetloborsky and Kamenushinsky massifs (Urals Platinum Belt)[END_REF]; solid-state deformation of the Galmoenan complex during ascent and emplacement, [START_REF] Burg | Translithospheric Mantle Diapirism: Geological Evidence and Numerical Modelling of the Kondyor Zoned Ultramafic Complex (Russian Far-East)[END_REF].

In addition, syn-emplacement deformation is compatible with a translithospheric mantle diapirism mechanism [START_REF] Bodinier | The Konder « Anorogenic » Ultramafic Massif (Aldan Shield) : Impact of a Translithospheric Mantle Diapir at the Earth's Surface ? Fourth International Workshop on Orogenic Lherzolite and Mantle Processes[END_REF][START_REF] Burg | Translithospheric Mantle Diapirism: Geological Evidence and Numerical Modelling of the Kondyor Zoned Ultramafic Complex (Russian Far-East)[END_REF]. conduit created within the dunite unit [START_REF] Augé | Primary platinum mineralization in the Nizhny Tagil and Kachkanar ultramafic complexes, Urals, Russia: a genetic model for PGE concentration in chromite-rich zones[END_REF]. In this model, the cavities created by physical processes were permanently fed by the magma going through the dunite body, the chromite being accumulated by convective processes (Augé et al, 2005, Fig. 14). Figure 10 illustrates this cavity model in the case of a Y-shaped Alaskan-type complex, with some field examples of different surface signatures of Alaskan-type complexes. It is here suggested that both high strain rate zones and the inclined edges of dunite cores of Alaskantype complexes favour at a smaller scale magma pathways and thus mineralizing processes within cavities.

Discussion and conclusion

Alaskan-type bodies have been classified according to their geometrical features at the surface and to their supposed structure at depth. Clustering of several of such bodies within typical distance of tens to hundreds of kilometres represents a specific observation which was attempted to be reproduced through numerical experiments. Based on the emplacement mechanism defined by the ascent of hot magma from the base of the lithosphere through a dyke-like vertical conduit, numerical models have been performed to test the effect of different parameters (thermal regimes, conduit diameter, and rheological parameters) on the mechanical response of the upper crust. After a series of one hundred high resolution experiments, it appears that the only way to keep a vertical dyke-like conduit, while weakness zones develop in the upper crust, is to consider a strong lower crustal rheology and some particular softening mechanisms, including tensile failure and strong strain dependence, in the brittle and ductile domains.

A C C E P T E D M A N U S C R I P T

Using ductile softening rules controlled by minimum and maximum accumulated strain values, a range of "Alaskan-type mechanical responses" have been obtained: (i)

cluster of cylindrical and funnel-shaped structures, including laterally shifted bodies, separated by 10-40 km, with topographic uplift (Figure 6); (ii) cluster of cylindrical structures, separated by 10-30 km, possibly connected at depth, with pronounced surface subsidence (Figure 7); (iii) cluster of Y-shape structures, separated by 25-55 km, with no significant erosion/subsidence (Figure 7a), and (iv) cluster of funnel-shaped structures, separated by 10-32 km, with pronounced surface subsidence (Figure 7b).

Because of the restricted 2D numerical procedure, a detailed comparison of geometrical aspects of the clustering process is not feasible. In addition, the 2D Cartesian geometry involves a "dyke-like" conduit rather than a "pipe-like" conduit. Consequently, one could simply observe from our simulations that several separated weakness zones can develop within the upper crust over a distance of 100-150 km. The approach used in this study corresponds to a series of large-scale thermomechanical numerical models applied to metallogenic issues. It has been demonstrated by several authors that the investigation of magma emplacement mechanisms, and more generally of large-scale geodynamic processes, could bring new insights into the understanding of mineralization processes (e.g. [START_REF] Barley | The late Archaean bonanza: metallogenic and environmental consequences of the interaction between mantle plumes, lithospheric tectonics and global cyclicity[END_REF][START_REF] De Boorder | Late Cenozoic mineralization, orogenic collapse and slab detachment in the European Alpine Belt[END_REF][START_REF] Guillou-Frottier | Genetic links between ash-flow calderas and associated ore deposits as revealed by large-scale thermo-mechanical modeling[END_REF][START_REF] Pirajno | Hotspots and mantle plumes: global intraplate tectonics, magmatism and ore deposits[END_REF]Bierlin et al., 2006;[START_REF] Burov | Plume headlithosphere interactions near intra-continental plate boundaries[END_REF][START_REF] Tosdal | From source to sinks in auriferous magmatichydrothermal porphyry and epithermal deposits[END_REF][START_REF] Begg | Lithospheric, cratonic, and geodynamic setting of Ni-Cu-PGE sulphide deposits[END_REF][START_REF] Guillou-Frottier | Plume-induced dynamic instabilities near cratonic blocks: implications for P-T-t paths and metallogeny[END_REF][START_REF] Bertrand | Distribution of porphyry copper deposits along the western Tethyan and Andean subduction zones: Insights from a paleotectonic approach[END_REF]. Here, we have studied a peculiar magmatic process that involves thermal and mechanical processes at both large (mantle) and small (crust) scales. At these scales, rock properties spanned a wide range of temperatures and pressures, implying that particular caution had to be paid when rheological properties were assessed.

Since rheological parameters seem to control the way the upper crust is fractured, there is a priori no reason to assign a particular geodynamic context (e.g. subduction zone)

to Alaskan-type bodies. This is consistent with the distinct geodynamic contexts where they occur (mainly in intra-continental and subduction-related settings). However, water-rich ultramafic magmas are probably more easily stored at the base of a continental lithosphere adjacent to a subduction zone.

The fast ascent of ultramafic magma through the lithosphere, which has been invoked by many authors (e.g. [START_REF] Beck | Dunite channels as viable pathways for mare basalts transport in the deep lunar mantle[END_REF][START_REF] Chen | Geochemistry and Os-Nd-Sr isotopes of the Gaositai Alaskan-type ultramafic complex from the northern North China craton: implications for mantle-crust interaction[END_REF] has been reproduced with a strong viscosity reduction due to strong strain softening. It can be underlined that our ascent velocity is ~50-100 times higher than in Gerya and Burg's models (2007), whereas effective viscosities remain similar. The role of strain softening in viscosity reduction may thus be critical in the dynamical features of ultramafic magmatism.

Among the several unresolved issues which were presented in the introduction, we did not tackle directly the problem dealing with the restricted age range of Alaskan-type complexes. However, it was shown that the favourable rheological conditions required a strong lower crust rather than a weak and ductile lower crust. The possibility of a strong lower crust would thus exist on Earth for only ~500 Myrs, suggesting that a weak ductile crust was predominant prior to this time. [START_REF] Connolly | Computation of phase equilibria by linear programming: a tool for geodynamic modelling and its application to subduction zone decarbonation[END_REF]. Right: strain rate field in case of thermodynamically consistent density (at 30 kyr). At this stage, the magmatic channel has been "destroyed" (laterally eroded). The strain rate field illustrates a nearly vertical conduit in the upper crust (upper cut-off is 5 10 -12 sec -1 , real maximal strain rates in the model are on the order of 10 -10 sec -1 ).

A C C E P T E D M A N U S C R I P T [START_REF] Nazimova | The Galmoenan massif, Kamchatka, Russia: geology, PGE mineralization, applied mineralogy and beneficiation[END_REF], Goodnews Bay [START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF], Kachkanar [START_REF] Augé | Primary platinum mineralization in the Nizhny Tagil and Kachkanar ultramafic complexes, Urals, Russia: a genetic model for PGE concentration in chromite-rich zones[END_REF], Konder [START_REF] Malitch | Assessment of the platinum potential of clinopyroxenite-dunite massifs[END_REF], Nizhny Tagil [START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF][START_REF] Augé | Primary platinum mineralization in the Nizhny Tagil and Kachkanar ultramafic complexes, Urals, Russia: a genetic model for PGE concentration in chromite-rich zones[END_REF],

Polaris [START_REF] Nixon | Geology and platinum-groupelement mineralization of Alaskan-type ultramafic complexes in British Columbia[END_REF], Tulameen [START_REF] Nixon | Geology and platinum-groupelement mineralization of Alaskan-type ultramafic complexes in British Columbia[END_REF][START_REF] Johan | The geology, geochemistry, mineralogy and mineral beneficiation of platinum-group elements[END_REF], Uktus [START_REF] Garuti | Chromite composition and platinum-group mineral assemblage in the Uktus Uralian-Alaskan-type complex (Central Urals, Russia)[END_REF], Wrede [START_REF] Nixon | Geology and platinum-groupelement mineralization of Alaskan-type ultramafic complexes in British Columbia[END_REF], Guli [START_REF] Malitch | Assessment of the platinum potential of clinopyroxenite-dunite massifs[END_REF], Kytlym (Zacarini et al., 2011).

Normalization values are taken from [START_REF] Barnes | The use of mantle normalization and metal ratios in discriminating between the effects of partial melting, crystal fractionation and sulphide segregation on platinum-group elements, gold, nickel and copper: examples from Norway[END_REF]

  and 3c) we have made a conventional assumption on the strength profile of the background lithosphere assuming initial geotherm yielding 60-70 mW/m² at surface (thermotectonic age of 100 Ma), granite-dominated crust, quartz-diorite lower crust and dry olivine mantle, and no ductile softening. We hence use olivine rheological parameters for mantle and quartz-granite rheology for the crust (see Supplementary material). The upper boundary condition is free surface, the bottom condition corresponds to Winkler basement (hydrostatic), and the lateral conditions are zero horizontal velocities. The thermal boundary conditions are zero temperature at surface, 1400°C at the bottom of the model (300 km depth), and zero lateral outflows. The phase changes are computed for the following compositions: mantleperidotite; crustgranite; magmapartially molten peridotite.

For

  these reasons, further series of high resolution experiments tempted to find the rheological conditions under which the magmatic material can mechanically reach the uppermost crustal layers without destroying the integrity of the magmatic conduit. The major parameters that have been tested within these experiments were: (1) the lower crustal rheology; (2) mantle rheology and geotherm; (3) softening mechanisms in the brittle and ductile parts; (4) width of the magmatic conduit; (5) magma temperature.

  Griffith failure can be controlled by varying T 0 (T 0 = 0 means permanent failure, T 0 = 1000 MPa means practically no tensile fracturing). As mentioned, this criterion applies when pressure becomes positive, in case of dilatation, basically in the uppermost rock layers (first kilometers in depth) where tensile stresses created by magma upraise and upward bending of the uppermost layers can negate lithostatic compressive pressure due to rock weight. Created tensile faults are "mode I" faults, i.e. these faults are normal to the tension axis while the Mohr-Coulomb faults are "mode II" faults. In deeper rocks total pressure is practically A C C E P T E D M A N U S C R I P T always negative due to the rapidly growing weight of the rock column (about 25-30 MPa per km), so Mohr-Coulomb failure prevails. The experiments detailed below show representative cases in which magmatic material arrives to break through the brittle crust and find its way to the surface. The series of experiments demonstrate the effects of different strain softening mechanisms in the ductile part. For the brittle part, we use linear cohesion softening, where cohesion drops from 20 MPa to 0 in the strain interval between 0.1 and 0.5. The used softening laws account for cumulative effects of grain size reduction, variations in fluid content in the deforming rock, partial melting and other effects mentioned in the previous sections. It is noteworthy that the choice of the form of the softening rule basically has no quantifiable experimental justification due to multiple processes involved. For example, Huismans and Beaumont (2003) assumed

  At a smaller scale, a detailed study of the Nizhny Tagil Complex shows that the Platinum Group Minerals (PGM) are systematically included in (or attached to) chromite crystals and are dominated by Ir alloys and Pt alloys, which explains the shape of the A C C E P T E D M A N U S C R I P T spidergram in Figure 9a. There is a correlation between the PGE content of a given sample and its chromite content. A model has been proposed to explain this relationship, where early crystalized PGM (platinum-group minerals) will be trapped by chromite. The subsequent concentration and accumulation of chromite crystals will cause the PGE enrichment in chromitites. Study of morphology of chromite concentration in the Nizhny Tagil Complex permitted to propose a genetic model where chromite accumulates in cavities along magma

  Besides reproducing individual or clustered weakness zones(Figures 4, 5 and 7), the results shown in Figure8raised an interesting question dealing with the apparent relationship between locations of high strain rate values and mineralization features reported in Figure9. Indeed, supposed Y-shaped geological structuresor twin Alaskan-type bodies(Galmoenan, Goodnews Bay, Kachkanar, etc) show the highest PGE concentrations, and Yshaped numerical bodies have the highest strain rate values. On the opposite, supposed cylindrical geological structuresor single Alaskan-type bodies (Konder, Nizhny Tagil, Tulameen, etc) have the lowest PGE concentrations, while cylindrical numerical bodies exhibit the lowest strain rate values.Obviously, the lack of data on Alaskan-type intrusion prevents to go further in the interpretation of this relationship. Nevertheless, scarce examples of syntectonic funnelshaped mafic/ultramafic intrusions exhibit similar geometrical relationships between both the location of melt injections and ore location (NiCu ± PGE mineralization), and the high strain rate zones. For example, the Variscan Aguablanca gabbronorite intrusion (SW Iberia,[START_REF] Romeo | 3D gravity modelling of the Aguablanca Stock, tectonic control and emplacement of a Variscan gabbronorite bearing a Ni-Cu-PGE ore, SW Iberia[END_REF]) is a funnel-shaped mafic intrusion emplaced in a major shear zone. The Ni-Cu-A C C E P T E D M A N U S C R I P TPGE mineralization is hosted by a funnel-like magmatic breccia located in the upper part of the funnel, close to the northern border of the intrusion. The mineralization is related to a final magmatic pulse injected within the shear zone (i.e. in the highest strain rate zone) controlling the mafic pluton emplacement. As a second example, the Permian Huangshandong ultramafic/mafic syntectonic intrusion (SW China,[START_REF] Branquet | Synkinematic mafic/ultramafic sheeted intrusions: Emplacement mechanism and strain restoration of the Permian Huangshan Ni-Cu ore belt (Eastern Tianshan, NW China)[END_REF] is also a funnelshaped intrusion with a sigmoid map shape emplaced within a major shear zone. NiCu mineralization is located close to the border and to the bottom of the intrusion where several coeval melt injections (magmatic mingling) are located and controlled by the shear zone bounding the syntectonic intrusion. These analogies also suggest a control of high strain rate zones on successive magmatic injections and mineralization location.

Figure 1 :

 1 Figure captions

Figure 2 :

 2 Figure 2: Examples of 20 Alaskan-type complexes: a) single circular or elliptical morphology;

Figure 3 :

 3 Figure 3: Numerical models of diapiric ascent of hot magma from the base of the lithosphere

Figure 4 :

 4 Figure 4: Case 1 ( 0 = 0.2,  1 = 0.5 ,  = 0.99900). Strong perturbations in the upper crust

Figure 5 :Figure 7 :

 57 Figure 5: Case 2 ( 0 = 0.1,  1 = 1.0 ,  = 0.99900). Strain rate pattern, showing regularly

Figure 8 :

 8 Figure 8: Summary of the individual geometries of the obtained crustal conduits: weakness

Figure 9 :

 9 Figure 9: PGE concentrations measured on chromitites, from the following Alaskan-type

  . a) typical signatures of Alaskantype complexes, defining a "M pattern" (grey area), with highest concentrations in Ir and Pt; b) data from anomalously mineralized chromites (see text); c) data from single Alaskan-type bodies, presenting the lowest concentrations; d) data from twin bodies, gathering the highest values.

Figure 10 :

 10 Figure 10: Sketch of an Alaskan-type complex with a Y-shape geometry (left), as inferred

Figure 1 A

 1 Figure 1

Figure 4 A

 4 Figure 4