Gaëtan Launay

Yannick Branquet

Stanislas Sizaret

Laurent Guillou- Frottier

Eric Gloaguen

Laurent Guillou-Frottier

How greisenization could trigger the formation of large vein-and-greisen Sn-W deposits: a numerical investigation applied to the Panasqueira deposit

Keywords: Panasqueira W-Sn deposit, Greisenization, reaction-enhanced permeability, fluid flux, numerical modeling

The formation of large tin-tungsten (Sn-W) deposits around granitic intrusions requires the circulation of large volumes of fluids within permeable structures. Half of the world's tungsten production originates from highly mineralized veins above granitic intrusions and from the altered part of the granite (the greisen), whose formation results from intense fluid-rock interactions. During greisenization processes, mineral reactions involve a decrease in the rock volume and thus an increase in porosity and permeability. To understand the complex fluid-rock interactions leading to the formation of large Sn-W ore deposits, we conducted numerical modeling accounting for magmatic fluid production and realistic permeability changes due to granite alteration and overpressure in the hosting rocks. The water/rock ratio is computed to constrain the rate of greisenization and therefore the porosity and permeability evolution laws. Four model results are presented: with and without fluid production exsolved from the granitic magma, and with and without dynamic reaction-enhanced permeability. The formation of greisen is reproduced, and greisen thickness reaches 200 m for the more sophisticated model. The interplay between greisenization and fluid production creates zones of overpressure above the granite that could localize the permeable structures such as the veins swarm observed at Panasqueira. Dynamic permeability promotes high fluid velocity and intense fluid-rock exchanges that could result in the formation of large ore deposits by enhancing mass transfer within and above granitic intrusions.

Introduction

The vein, stockwork, and greisen Sn-W deposits form parts of magmatic-hydrothermal systems that provide an important part of the world's tungsten production [START_REF] Werner | International strategic mineral issues summary report-Tungsten[END_REF] and represent an important source of tin. These deposits generally include: (i) disseminated mineralization within a massive greisen characterized by low grade and high tonnage and (ii) mineralized quartz veins usually characterized by high grade and low tonnage [START_REF] Taylor | Geology of tin deposits[END_REF][START_REF] Taylor | Pervasive hydrothermal alteration in tin-bearing granite and implications for the evolution of ore-bearing magmatic fluids[END_REF][START_REF] Mlynarczyk | San Rafael, Peru: Geology and structure of the world's richest tin lode[END_REF]. The formation of veins is strongly controlled by the regional structures and the fluid pressure conditions that promote the opening of the permeable structures [START_REF] Sibson | Earthquake rupturing as a mineralizing agent in hydrothermal systems[END_REF][START_REF] Sibson | Preparation zones for large crustal earthquakes consequent on fault-valve action[END_REF][START_REF] Cox | Coupling between deformation, fluid pressures and fluid flow in ore-producing hydrothermal environments[END_REF]2010). Massive greisens are generally present in the apical portions of granite intrusions and result from strong interactions between granitic rocks and hydrothermal fluids during the magmatic-hydrothermal transition [START_REF] Stemprok | Greisenization (a review)[END_REF][START_REF] Pollard | Genetic Modelling of Greisen-Style Tin Systems[END_REF][START_REF] Pirajno | Greisen systems[END_REF][START_REF] Černý | Granite-Related Ore Deposits[END_REF]. The formation of vein and greisen deposits implies the transfer and the focusing of a large volume of mineralizing fluids [START_REF] Dolejs | Quantitative characterization of hydrothermal systems and reconstruction of fluid fluxes: the significance of advection, disequilibrium, and dispersion[END_REF]. Moreover, a recent fluid inclusion study performed by [START_REF] Korges | Depressurization and boiling of a single magmatic fluid as a mechanism for tin-tungsten deposit formation[END_REF] in the vein-and-greisen system of Zinnwald (Erzgebirge, in between the Czech Republic and Germany) suggests that the formation of mineralized veins and greisen are contemporaneous and associates with the boiling of a single magmatic fluid.

The advection constitutes a privileged mode of heat and mass transport, the efficiency of which directly depends on rock permeability and the driving forces involved during fluid flow [START_REF] Phillips | Flow and reactions in permeable rocks[END_REF][START_REF] Hedenquist | The role of magma in the formation of hydrothermal ore deposits[END_REF][START_REF] Cox | Coupling between deformation, fluid pressures and fluid flow in ore-producing hydrothermal environments[END_REF]). According to the geochemical study performed by Cernushi et al. (2018), the mineralization stage around granitic intrusions would occur over a very short period (less than 35 kyr for large-scale economic porphyry copper deposits). Similarly, the formation of large-scale ore deposits around granitic intrusions would require a very efficient heat and mass transfer process.

Magmatic-hydrothermal mineralized systems were investigated in numerous studies involving numerical simulation describing and constraining the fluid flow pattern in and around granitic intrusions [START_REF] Norton | Thermal aspects of ore deposition[END_REF][START_REF] Cathles | How long can a hydrothermal system be sustained by a single intrusive event?[END_REF][START_REF] Gerdes | Convective fluid flow through heterogeneous country rocks during contact metamorphism[END_REF][START_REF] Cui | Controls of layered and transient permeability on fluid flow and thermal structure in contact metamorphic aureoles, with application to the Notch Peak aureole Utah[END_REF][START_REF] Eldursi | Numerical investigation of transient hydrothermal processes around intrusions: Heat-transfer and fluid-circulation controlled mineralization patterns[END_REF]2018). These models permit to constrain in a simplified way (i.e., with static permeability values) the time evolution of the fluid flow patterns and, in some cases, to delineate in time and space the most favorable zones for mineralization to occur. However, permeability is a dynamic parameter that is mainly affected by interplays between the stress state, the rock alteration and fluid pressure conditions [START_REF] Sibson | High-angle reverse faults, fluid-pressure cycling, and mesothermal gold-quartz deposits[END_REF][START_REF] Manning | Permeability of the continental crust: implications of geothermal data and metamorphic systems[END_REF][START_REF] Cox | The application of failure mode diagrams for exploring the roles of fluid pressure and stress states in controlling styles of fracture-controlled permeability enhancement in faults and shear zones[END_REF][START_REF] Ingebritsen | The physical hydrogeology for ore deposits[END_REF][START_REF] Ingebritsen | Crustal permeability: Introduction to the special issue[END_REF]. Effects of rock deformation and fracturing on permeability and hence on fluid flow were investigated in numerous structural and experimental studies [START_REF] Sibson | Earthquake rupturing as a mineralizing agent in hydrothermal systems[END_REF]1988;[START_REF] Zhang | The influence of room temperature deformation on porosity and permeability in calcite aggregates[END_REF]Tenthorey and Fitz Gerald, 2006;[START_REF] Coelho | Permeability of sheeted dykes beneath oceanic ridges: Strain experiments coupled with 3D numerical modeling of the Troodos Ophiolite, Cyprus. Tectonophysics[END_REF].

These studies have demonstrated that deformation transiently enhances the permeability of the hydrothermal systems by the formation of permeable structures like shear zone and fractures. Fluid overpressure conditions related to fluid production are another important parameter that promotes the opening of permeable structures enhancing the rock permeability by factors of 100 to 1000 [START_REF] Zhang | The influence of room temperature deformation on porosity and permeability in calcite aggregates[END_REF][START_REF] Ingebritsen | Crustal permeability: Introduction to the special issue[END_REF].

The hydrologic effects of these different feedbacks were investigated by numerical modeling in cases of geothermal reservoirs related to magmatic-hydrothermal ore deposits by [START_REF] Weis | Porphyry-copper ore shells form at stable pressure-temperature fronts within dynamic fluid plumes[END_REF] and [START_REF] Weis | The dynamic interplay between saline fluid flow and rock permeability in magmatic hydrothermal systems[END_REF].

These numerical studies have highlighted that the combination of dynamic permeability with the stress state and the magmatic fluid production permits to reproduce the geometry of ore shells commonly observed in copper porphyry deposits [START_REF] Weis | Porphyry-copper ore shells form at stable pressure-temperature fronts within dynamic fluid plumes[END_REF].

Depending on the spatial organization of fluid flow, the pressure and temperature conditions and the chemical characteristics of fluids, hydrothermal alteration is associated with the dissolution, precipitation, and mineralogical transformation of rocks that affect directly the permeability [START_REF] Putnis | Mineral replacement reactions: from macroscopic observations to microscopic mechanisms[END_REF]Jamtveit et al., 2009;[START_REF] Putnis | Fluid-induced processes: metasomatism and metamorphism[END_REF][START_REF] Pollok | Mineral replacement reactions in solid solution-aqueous solution systems: Volume changes, reactions paths and endpoints using the example of model salt systems[END_REF][START_REF] Booden | Quantifying Metasomatism in Epithermal Au-Ag Deposits: A Case Study from the Waitekauri Area, New Zealand[END_REF][START_REF] Putnis | Replacement processes in the earth's crust[END_REF][START_REF] Jonas | The role of grain boundaries and transient porosity in rocks as fluid pathways for reaction front propagation[END_REF][START_REF] Putnis | Transient Porosity Resulting from Fluid-Mineral Interaction and its Consequences[END_REF][START_REF] Scott | Permeability Changes Resulting from Quartz Precipitation and Dissolution around Upper Crustal Intrusions[END_REF]. According to the volume changes of the solid phase induced by replacement reactions related to hydrothermal alteration, these mineralogical transformations can lead to the self-sealing of the permeable pathways causing the inhibition of fluid flow, or conversely result in porosity and permeability increase that enhances fluid flow and promote mass and heat transport. The greisenization responsible for the massive greisen formation is a hydrothermal alteration leading to partial to total replacement of early magmatic feldspars and biotites of granite by a late magmatic-hydrothermal quartz muscovite assemblage. An experimental and petrophysical study performed by [START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF] has demonstrated that the replacement of feldspars by muscovite during the greisenization is accompanied by a volume decrease of the mineral phases of about 8.5%. This volume decrease will induce porosity generation and hence permeability increase that can reach 3 orders of magnitude.

Massive greisens are common in the upper part of granitic intrusions related to Sn-W mineralized systems, and their thickness can reach ~ 200-300 meters like at Cinovec, East Kemptville and Panasqueira [START_REF] Halter | Origin and evolution of the greisenizing fluid at the East Kemptville tin deposit[END_REF][START_REF] Jarchovsky | The nature and genesis of greisen stocks at Krásno, Slavkovský les area western Bohemia, Czech Republic[END_REF][START_REF] Laznicka | Giant metallic deposits. Future sources of industrial metals[END_REF]. The development of such a thick permeable pathway during the magmatic-hydrothermal activity responsible for the formation of Sn-W deposits can significantly enhance fluid flow and mass transport. However, the effects of this feedback between the greisenization and permeability on fluid flow in vein and greisen deposits need to be constrained. Based on the coupling of physical equations that describe the physical hydrology of fluid flow and heat and mass transport in hydrothermal systems, the numerical modeling approach represents a powerful tool to understand and quantify the controlling processes leading to the formation of hydrothermal ore deposits, as well as the development of hydrothermal alteration zones.

As mentioned above, the potential mechanisms susceptible to affect the permeability and fluid flow during the formation of vein and greisen deposits are: (i) the greisenization of granite that can enhance the permeability of granite roof, promote fluid flow and facilitate the expulsion of magmatic fluids released during the granite crystallization at depth and (ii) the opening of permeable structures in response to fluid overpressure conditions induced by the expulsion of magmatic fluids. In this study, we explore by numerical modeling the interplays between the permeability changes induced by the greisenization, the production of magmatic fluid during granite crystallization, and the opening of permeable structures induced by potential fluid overpressure conditions. The model geometry and parameters are based on geological and petrophysical observations (porosity and permeability measurements) from the world-class Panasqueira W-Sn-(Cu) ore deposit. Although Panasqueira is an atypical granite-related W-Sn system in terms of vein geometry (subhorizontal), it represents a reference system to study magmatic-hydrothermal processes leading to the formation of large vein and greisen deposits. We present results from 2D numerical modeling of the time evolution of fluid flow patterns and fluxes to evaluate the influence of permeability changes during greisenization processes.

The Panasqueira W-Sn-(Cu) deposit

The Panasqueira W-Sn-(Cu) deposit is located in the Central Iberian Zone (CIZ) in the middle part of the Iberian Variscan Belt (Fig. 1a) [START_REF] Julivert | Mapa Tectonico de la Peninsula Iberica y Baleares E[END_REF]. This tectonic domain is marked by occurrences of numerous granitic intrusions mainly emplaced into a thick monotonous Neoproterozoic schist-greywacke sequence during the latest stages of the Variscan orogeny [START_REF] Dias | U-Pb zircon and monazite geochronology of syn-to post-tectonic Hercynian granitoids from the central Iberian Zone (northern Portugal)[END_REF][START_REF] Castro | Palaeozoic Magmatism[END_REF]. The CIZ is particularly rich in Sn-W deposits and showings, whose formation is related to the magmatic-hydrothermal processes induced by the emplacement of a large volume of granitic bodies. These deposits comprise mostly vein and greisen systems like Regoufe, Gois, and Panasqueira, which were/are mainly exploited for the W mineralization.

The world-class Panasqueira W-Sn-(Cu) deposit belongs to the Beira Baixa province located in the south of the Serra da Estrela granitic Massif. This region is largely composed of a schist-greywacke sequence affected by tight and upright folds and by a low-grade greenschist metamorphism during the Variscan orogeny.

The presence of spotted schist delineates a thermal metamorphism aureole, which betrays the presence of a non-outcropping granitic intrusion at depth [START_REF] Clark | Preliminary study of the temperatures and confining pressures of granite emplacement and mineralization[END_REF][START_REF] Kelly | Geologic, fluid inclusion and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF][START_REF] Bussink | Geochemistry of the Panasqueira Tungsten-Tin Deposit[END_REF]. The surface distribution of spotted schist, the drill holes data over the Panasqueira district, and gravimetric survey show that the granite body is elongated in the NW-SE direction (7.5 x 4.5 km). Its thickness does not exceed 2 km (Fig. 1b) [START_REF] Hebblethwaite | A report on the study of dilation patterns within the Panasqueira ore body: Barroca Grande[END_REF][START_REF] Thadeu | Geologia do couto mineiro da Panasqueira[END_REF][START_REF] Kelly | Geologic, fluid inclusion and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF][START_REF] Clark | Preliminary study of the temperatures and confining pressures of granite emplacement and mineralization[END_REF].

A cupola intersected by the underground mining work and the upper part of this granite were altered into a massive greisen during fluid flow related to the incipient stage of the magmatic-hydrothermal activity at Panasqueira (Fig. 1b) [START_REF] Thadeu | Geologia do couto mineiro da Panasqueira[END_REF][START_REF] Clark | Preliminary study of the temperatures and confining pressures of granite emplacement and mineralization[END_REF][START_REF] Kelly | Geologic, fluid inclusion and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF][START_REF] Bussink | Geochemistry of the Panasqueira Tungsten-Tin Deposit[END_REF][START_REF] Launay | Genetic relationship between greisenization and Sn-W mineralization in vein and greisen deposits: Insights from the Panasqueira deposit (Portugal)[END_REF].

This massive greisen can be characterized by porous and microporous textures emphasizing the development of porosity during the replacement reactions related to the greisenization processes. Petrophysical measurements performed by [START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF] highlight: (i) a significant porosity generation (~8.5%) during the greisenization that can be also observed in other mineralized greisen systems and (ii) a significant permeability increase (more than 3 orders of magnitude) induced by this porosity generation. The occurrences of metals-bearing minerals like cassiterite and chalcopyrite in the pores generated during the greisenization suggest that the development of these permeable pathways could promote and enhance fluid flow in greisen systems during the mineralizing events.

The economic W-Sn mineralization consists of a dense swarm of sub-horizontal quartz veins mostly hosted by the metasedimentary host rock (Beira schist). The ore zone is spatially associated with the greisen cupola and extends over an area of 6 km² (Fig. 1b and1c). The presence of geochemical tungsten anomalies centered on the cupola [START_REF] Oosterom | Lithogeochemical studies of aureoles around the panasqueira tin-tungsten deposit[END_REF] suggests an important role of this cupola on fluid flow organization (Fig. 1b), as confirmed by [START_REF] Launay | Deciphering fluid flow at the magmatic-hydrothermal transition: A case study from the world-class Panasqueira W-Sn-(Cu) ore deposit (Portugal)[END_REF]2021). The Panasqueira district also has late NNW-SSE and ENE-WSW sub-vertical strike-slip fault systems crosscutting and postdating the mineralized vein systems of Panasqueira re-activated during the Alpine orogeny (Fig. 1b) [START_REF] Thadeu | Geologia do couto mineiro da Panasqueira[END_REF][START_REF] Kelly | Geologic, fluid inclusion and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF][START_REF] Foxford | Extensional failure and hydraulic valving at Minas da Panasqueira, Portugal: Evidence from vein spatial distributions, displacements and geometries[END_REF]. Consequently, these faults were not active during fluid flow related to the W-Sn mineralization and hence are not considered in the present work.

The formation depth of the Panasqueira deposit is abundantly debated. [START_REF] Cathelineau | High pressure and temperatures during the early stages of tungsten deposition at Panasqueira revealed by fluid inclusions in topaz[END_REF] reviewed P/T estimates for the genesis of this ore deposit. Based on new fluid inclusion data, these authors propose a formation depth for the Panasqueira deposit as deep as 8-10 km with a high geothermal gradient of about 60°C/km. However, several lines of evidence invalidate such a deep scenario: (i) the CIZ domain has been weakly thickened during the Variscan orogen and no allochthonous thrust sheet covered the Panasqueira area at that time. As a result, the only folding of Paleozoic sedimentary sequences cannot account for such deep subsidence to 8-10 km [START_REF] Martínez Poyatos | Imaging the crustal structure of the Central Iberian Zone (Variscan Belt): The ALCUDIA deep seismic reflection transect[END_REF][START_REF] Díez-Fernández | Extensional orogenic collapse captured by strike-slip tectonics: Constraints from structural geology and U-Pb geochronology of the Pinhel shear zone (Variscan orogen, Iberian Massif)[END_REF]; (ii) the Beira schist hosting the ore deposits suffered a very weak greenschist facies metamorphism (e.g., fossils and sedimentary structures are preserved and regional chlorite isograde is not reached far from the thermal aureole [START_REF] Clark | Preliminary study of the temperatures and confining pressures of granite emplacement and mineralization[END_REF][START_REF] Kelly | Geologic, fluid inclusion and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF][START_REF] Bussink | Geochemistry of the Panasqueira Tungsten-Tin Deposit[END_REF]. Such a low metamorphic grade is not compatible with subsidence to 8-10 km; (iii) at a depth of 8-10 km, permeabilities are very low and promote dominant diffusion vs. advection for mass and heat transfers (e.g., [START_REF] Manning | Permeability of the continental crust: implications of geothermal data and metamorphic systems[END_REF], This leads [START_REF] Cathelineau | High pressure and temperatures during the early stages of tungsten deposition at Panasqueira revealed by fluid inclusions in topaz[END_REF] and [START_REF] Carocci | Incipient Wolframite Deposition at Panasqueira (Portugal): W-Rich Rutile and Tourmaline Compositions as Proxies for the Early Fluid Composition[END_REF] to propose the involvement of deep-rooted high-permeability shear zones in the Panasqueira area. Such shear zones have not been mapped nor observed in the area; (iv) using different approaches (fluid inclusions, structural analysis, and rock mechanics) many researchers working on the Panasqueira deposit proposed that the granite body emplaced at a shallow crustal level, at about 3 km (e.g., [START_REF] Kelly | Geologic, fluid inclusion and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF][START_REF] Bussink | Geochemistry of the Panasqueira Tungsten-Tin Deposit[END_REF][START_REF] Polya | Chemistry of the main-stage ore-forming fluids of the Panasqueira W-Cu-(Ag)-Sn deposit, Portugal: implications for models of ore genesis[END_REF][START_REF] Foxford | Extensional failure and hydraulic valving at Minas da Panasqueira, Portugal: Evidence from vein spatial distributions, displacements and geometries[END_REF]. The neighboring Variscan Argemela rare-metal granite (13 km east of Panasqueira) displays a subvolcanic texture (snow-ball textured quartz) also indicating a shallow emplacement level (1-2 km, [START_REF] Michaud | Magmatic fractionation and the magmatic-hydrothermal transition in rare metal granites: Evidence from Argemela (Central Portugal)[END_REF]. Based on those critical points and following our former published works on the Panasqueira deposits system [START_REF] Launay | Deciphering fluid flow at the magmatic-hydrothermal transition: A case study from the world-class Panasqueira W-Sn-(Cu) ore deposit (Portugal)[END_REF][START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF], we favor a granite roof emplacement depth of about 3 km.

The cumulative vein thickness along the vertical sections of drill holes emphasizes that mineralized veins of Panasqueira form a regular band of about 100-250 m of thickness above the granite roof (Fig. 1d, after [START_REF] Foxford | Extensional failure and hydraulic valving at Minas da Panasqueira, Portugal: Evidence from vein spatial distributions, displacements and geometries[END_REF]. This distribution suggests that the fluid responsible for the mineralization of Panasqueira has preferentially flowed through a regular horizontal permeable fractured zone. From this observation, [START_REF] Foxford | Extensional failure and hydraulic valving at Minas da Panasqueira, Portugal: Evidence from vein spatial distributions, displacements and geometries[END_REF] proposed a structural model of vein formation involving extensional failure induced by supralithostatic fluid pressure in a compressive crustal regime (σ 3 = σ v ). Recent results based on the growth bands study of tourmaline are consistent with the implication of fluid overpressure during the incipient stages of the vein opening [START_REF] Launay | Deciphering fluid flow at the magmatic-hydrothermal transition: A case study from the world-class Panasqueira W-Sn-(Cu) ore deposit (Portugal)[END_REF]. However, the expulsion of magmatic fluids at lithostatic pressure conditions during the granite crystallization can constitute an important mechanism during the opening of the vein swarm and need to be explored by numerical modeling. The formation of the dense vein system was related to the circulation of a large volume of fluids estimated at ~1000 km 3 [START_REF] Polya | Chemistry of the main-stage ore-forming fluids of the Panasqueira W-Cu-(Ag)-Sn deposit, Portugal: implications for models of ore genesis[END_REF] for the W-Sn mineralization stages. This implies efficient processes of fluid flow and hence sufficient permeable structures and lithological facies to permit the transfer of this volume of fluid.

Model set-up and time-dependent properties

To investigate the role of greisenization on the formation of W-Sn deposits, we have performed a timedependent 2D numerical modeling of fluid flow and heat transfer around a granitic intrusion, integrating: (i) the depth-dependent permeability profile, (ii) the effect of fluid overpressure conditions on the permeability of the host rock in response to magmatic fluid expulsion, and (iii) effects of feedback between greisenization, porosity, and permeability on fluid flow. The acid fluids responsible for greisenization are mainly sourced from the crystallizing magmas implying that magmatic fluid production has to be accounted for in experiments. However, the convection of crustal-derived and/or meteoric fluids within host rocks are known to play an important role during Sn ore deposition associated with greisens near the granite roof (e.g., [START_REF] Fekete | Contrasting hydrological processes of meteoric water incursion during magmatic-hydrothermal ore deposition: An oxygen isotope study by ion microprobe[END_REF]. Therefore, assuming a sufficient acidity of fluids in host-rock to trigger greisenization, it is also worth testing configurations without magmatic fluid production. Thus, we have generated four different numerical models, involving two different types of magmatic-hydrothermal systems: (i) one with a "passive" granitic intrusion without fluid production and where fluid flow is only driven by convective processes, and

(ii) another one involving magmatic fluid production during granite crystallization and where overpressure fluid conditions constitute an additional driving force.

These models were performed with the finite element code Comsol Multiphysics™. This numerical code has been tested and validated by several benchmarks in cases of fluid flow around hot granite intrusions [START_REF] Eldursi | Numerical investigation of transient hydrothermal processes around intrusions: Heat-transfer and fluid-circulation controlled mineralization patterns[END_REF] and geothermal systems [START_REF] Garibaldi | Thermal anomalies and geological structures in the Provence basin: Implications for hydrothermal circulations at depth[END_REF][START_REF] Taillefer | Fault-Related Controls on Upward Hydrothermal Flow: An Integrated Geological Study of the Têt Fault System, Eastern Pyrénées (France)[END_REF][START_REF] Guillou-Frottier | On the morphology and amplitude of 2D and 3D thermal anomalies induced by buoyancy-driven flow within and around fault zones[END_REF] for which results were compared with those obtained from well-established codes.

Governing equations

In magmatic-hydrothermal systems, the main driving forces involved during fluid flow are: (i) the buoyancy force induced by variation of fluid density due to the change in temperature, pressure, and salinity conditions, and (ii) the production of magmatic fluids that affect the pressure gradient and hence the mass transport by advection. As the contribution of surface-derived fluids during ore deposition has not been established at Panasqueira, the topography-driven pressure gradients are not considered here. In this study, we consider several physical equations that were usually applied and coupled to describe the behavior of fluid flow in hydrothermal systems. The mass conservation equation is considered for a temperature and slightly pressure-dependent (i.e., near incompressible) fluid density within a saturated porous medium:

(1)

∂(𝜑.𝜌 𝑓 ) ∂𝑡 = -∇.(𝜌 𝑓 .𝑢) + 𝑄 𝐻 2 𝑂
where φ is the rock porosity, ρ f the fluid density (kg.m -3 ), t the time (s), u the fluid velocity (m.s -1 ), and Q H2O a fluid source term (kg.m -3 .s -1 ) in models corresponding to the magmatic fluid production.

Darcy's law was applied to calculate the fluid velocity according to the following equation:

(

2) 𝑢 = - 𝑘 𝜇 𝑓 (∇𝑃 -𝜌 𝑓 𝑔)
with k the rock permeability (m²), μ f the fluid dynamic viscosity (Pa.s), P the fluid pressure (Pa), and g the gravitational acceleration (m.s -2 ).

Heat transport is achieved by conduction and advection in a porous medium and is described for a single fluid phase by:

(3)

𝑄 𝐻 = 𝐶 𝑒𝑞( ∂𝑇 ∂𝑡 ) + ∇.( -𝜆 𝑒𝑞 ∇𝑇) + 𝜌 𝑓 𝐶𝑝 𝑓 .𝑢.∇𝑇
With Q H (W.m -3 ) the radiogenic heat source produced by the granite intrusion, T (K) the temperature, ρ f (kg.m - 3 ) the fluid density, Cp f (J.kg -1 .K -1 ) the specific heat capacity of fluid and C eq and λ eq are the weighted average volumetric heat capacity and equivalent thermal conductivity, defined in saturated porous media as:

(4) 𝐶 𝑒𝑞 = 𝜑.𝜌 𝑓 .𝐶𝑝 𝑓 + (1 -𝜑).𝜌 𝑟 .𝐶𝑝 𝑟 and:

(5) 𝜆 𝑒𝑞 = 𝜑.𝜆 𝑓 + (1 -𝜑).𝜆 𝑟 With ρ r (kg.m -3 ) the rock density, Cp r (J.kg -1 .K -1 ) the heat capacity of rocks, and λ r and λ f (W.m -1 .K -1 ) are, respectively, the thermal conductivities of the rocks and fluid. The latent heat, produced during the crystallization of the granite domain, is accounted by applying a doubled value of heat capacity (Cp r ) for the granite domain for temperatures above the solidus temperature (T > 600°C) [START_REF] Hayba | Multiphase ground water flow near cooling plutons[END_REF].

Model geometry, boundary conditions, and rock properties

The model geometry is based on geological and borehole data acquired at Panasqueira [START_REF] Hebblethwaite | A report on the study of dilation patterns within the Panasqueira ore body: Barroca Grande[END_REF][START_REF] Thadeu | Geologia do couto mineiro da Panasqueira[END_REF][START_REF] Kelly | Geologic, fluid inclusion and stable isotope studies of the tin-tungsten deposits of Panasqueira, Portugal[END_REF][START_REF] Wheeler | Technical report on the mineral resources and reserves of the Panasqueira mine, Portugal[END_REF][START_REF] Launay | Deciphering fluid flow at the magmatic-hydrothermal transition: A case study from the world-class Panasqueira W-Sn-(Cu) ore deposit (Portugal)[END_REF]2021). The chosen 2D geometry is inferred from the A-A' cross-section shown in Figures 1c and2. We assume a granite intrusion of a laccolith shape with a 6 km lateral extent and a height of 2 km. This granitic intrusion is asymmetric at 3km depth with a cupola in the granite roof at 2.7 km depth at the time of the intrusion. The granite intrusion emplaced in a metasedimentary host-rock domain of 20 km in width for 10 km in height. We assume an initial granite temperature of 850°C and an initial geothermal gradient of 30°C per kilometer in the host rock domain.

Initial hydrostatic fluid pressure is imposed throughout all the domains of the model. This situation is disputable compared to the natural conditions within the granite. During the early stages of granite emplacement, the magma is poorly crystallized, and the pores are weakly connected. Such conditions are favorable to lithostatic pressure regime within the granite body [START_REF] Chi | Hydrodynamic links between shallow and deep mineralization systems and implications for deep mineral exploration[END_REF][START_REF] Fournier | Hydrothermal processes related to movement of fluid from plastic into brittle rocks in the magmatic-epithermal environment[END_REF]. However, at the interface between the granite and the host rock, this colder and more crystallized zone is characterized by a transient pressure regime between lithostatic and hydrostatic. To avoid producing an initial abnormally high-pressure gradient at the interface surrounding the granite, we then assume an initial hydrostatic pressure regime throughout the model. This numerical choice represents a useful limiting case for discussing the influence of fluid production by granite on fluid pressure regime apart from all other influences. Moreover, the low permeability of the granite limits the fluxes within the volume of the pluton itself. The physical properties of the granite and the host rock domains were chosen according to the physical properties of granitic rocks and metasedimentary rocks (schist) commonly described in the literature (Table 1).

We assume a flat surface topography at the upper boundary with a fixed pressure of 1.015 × 10 5 Pa and a fixed temperature of 20 °C (Fig. 2). A temperature of 320 °C and impermeable conditions are fixed at the bottom contact. The lateral boundaries are assumed impermeable and thermally insulated (Fig. 2). The mesh comprises 24646 triangular elements with resolution ranging between element side lengths of 50 m to 200 m.

The elements were refined in the granite domain, along the upper boundary, and along the contacts between the granite domain and the host rock domain that constitute the most important part of the model.

Fluid properties

For the fluid phase, we assume a porous medium saturated by a single fluid phase composed of pure water. The fluid density was calculated at each time step as a function of the pressure and the temperature using an interpolated function based on the IAPWS-IFP97 equation of state for pure water defined by the International Association for the Properties of Water and Steam (IAPWS, [START_REF] Wagner | The IAPWS industrial formulation 1997forth thermodynamic properties of water and steam[END_REF]. From this database, an interpolated function (polynomial with 15 coefficients, as detailed in [START_REF] Guillou-Frottier | On the morphology and amplitude of 2D and 3D thermal anomalies induced by buoyancy-driven flow within and around fault zones[END_REF] was implemented as a parameter in the solver. Based on the analytical approximation of [START_REF] Kestin | Effect of the pressure on the viscosity of aqueous NaCl solutions in the temperature range 20-150 C[END_REF] and [START_REF] Rabinowicz | Two and three-dimensional modelling of hydrothermal convection in the sedimented Middle Valley segment, Juan de Fuca Ridge[END_REF], the dynamic viscosity evolution with temperature is given by: (6) 𝜇 𝑓 = 2.414 10 -5 .10 ( 247.8 𝑇 -140 )

with T the fluid temperature in K. As discussed by [START_REF] Clauser | Geothermal energy[END_REF], fluid viscosity is mainly controlled by temperature conditions and the effects of pressure can be neglected.

Magmatic fluid production

At each time step, the rate of magmatic fluid production is calculated as being proportional to the crystallization rate of the granitic intrusion between the liquidus (T l = 700°C) and the solidus (T s =600°C)

temperatures. For each element composing the granite domain characterized by temperatures ranging between the solidus and the liquidus, the volume of crystallized granite (V cryst ) was calculated from the melt fraction (F) according to the formulation of Caricchi and Blundy ( 2015):

(7) 𝑉 𝑐𝑟𝑦𝑠𝑡 (𝑡) = (1 -𝐹).𝐴.𝑑𝑦 with: (8) 𝐹 = 1 - 58.21 + 8.908sinh ( -4.831[ 𝑇 -𝑇𝑠 𝑇𝑙 -𝑇𝑠 ] + 2.257 ) 100
where V crys t (m 3 ) is the volume of crystallized granite at the considering time step, F the melt fraction, A (m²)

is the surface of the considering element and dy (m) the extension of the element in the third dimension; here we assume dy = 1 m to account for only the mass of fluid produced on the 2D section considered in our modeling. According to the compilation of water concentration data acquired from melt inclusions analyzed in minerals from worldwide granites [START_REF] Thomas | Water in granite and pegmatite-forming melts[END_REF], we assume that 5 wt% of fluid (pure water) was produced from the volume of crystallized granite (V crys t). The effective mass of fluid (M W (t))

produced at the considering time step was calculated by subtracting the mass of fluid produced during the previous time step:

(9) 𝑀 𝑊 (𝑡) = 0.05.𝜌 𝑟 .𝑉 𝑐𝑟𝑦𝑠𝑡 (𝑡) -𝑀 𝑊 (𝑡 -1)
with M W (t) (kg) the mass of fluid produced at the considering time step, ρ r (kg.m -3 ) the density of granite, V cryst (t) the volume of crystallized granite at the considering time step and M W (t-1) (kg) the mass of fluid produced during the previous time step. From this mass of produced fluid, the rate of fluid production is calculated according to the following expression:

(10) 𝑄 𝐻2𝑂 (𝑡) = 𝑀 𝑤 (𝑡). 1 𝑉.𝑑𝑡
where Q H2O (kg.m -3 .s -1 ) is the rate of fluid production at the considering time step, M W (t) (kg) is the mass of fluid produced by the considering element at the considering time step, V (m 3 ) is the volume of the considered element and dt (s) is the time difference between successive time steps. The mass conservation equation ( 1) describes the conservation of mass in a cell, i.e., the quantity of matter which enters in a cell is equal to that which comes out minus the accumulation term depending on porosity and fluid density variations through times (close to null here because the fluid is near incompressible). When the source term Q H20 is non-zero (magmatic fluid production herein), equation ( 1) is equilibrated with an increase of fluid velocity and therefore an increase of pressure gradient in the Darcy's equation (2).

Permeability of the host-rock domain

3.5.1

The background permeability of host-rocks

Compilation of permeability values estimated in tectonically active continental crust permitted to define a depth-dependent profile of permeability [START_REF] Manning | Permeability of the continental crust: implications of geothermal data and metamorphic systems[END_REF]. Due to the impermeable nature of the metasedimentary host rocks present at Panasqueira [START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF], we have defined, for depths exceeding 0.5 km, a function of permeability evolution modified from this depth-dependent function:

(11) log (𝑘) = -16 -3.2log (𝑧)
where k (m²) is permeability and z (km) is depth. Above 0.5 km, permeability was fixed at 10 -15 m² (Fig. 3a).

This function was applied in our models to define the background permeability values in the host rock domain for hydrostatic fluid pressure conditions.

Feedback between fluid overpressure conditions and permeability of host-rocks

Fluid overpressure conditions (P f > P hydrostatic ) constitute a key process to increase temporarily the rock permeability [START_REF] Sibson | High-angle reverse faults, fluid-pressure cycling, and mesothermal gold-quartz deposits[END_REF][START_REF] Cox | The application of failure mode diagrams for exploring the roles of fluid pressure and stress states in controlling styles of fracture-controlled permeability enhancement in faults and shear zones[END_REF][START_REF] Weis | The dynamic interplay between saline fluid flow and rock permeability in magmatic hydrothermal systems[END_REF]. Indeed, numerous experimental and field studies have emphasized that permeability can increase by 2 orders of magnitude under overpressure fluid conditions [START_REF] Faulkner | The effect of temperature, the nature of the pore fluid, and subyield differential stress on the permeability of phyllosilicate-rich fault gouge[END_REF][START_REF] Cox | Coupling between deformation, fluid pressures and fluid flow in ore-producing hydrothermal environments[END_REF][START_REF] Ingebritsen | Permeability of the continental crust: dynamic variations inferred from seismicity and metamorphism[END_REF][START_REF] Howald | Evidence for long-timescale (>103 years) changes in hydrothermal activity induced by seismic events[END_REF][START_REF] Ingebritsen | Crustal permeability: Introduction to the special issue[END_REF]. This feedback is mainly involved in hydraulic fracturing processes as fault-valve systems but can also be applied to explain the formation of veins in response to the expulsion of magmatic fluids in magmatic-hydrothermal deposits [START_REF] Weis | The dynamic interplay between saline fluid flow and rock permeability in magmatic hydrothermal systems[END_REF]. At Panasqueira, [START_REF] Foxford | Extensional failure and hydraulic valving at Minas da Panasqueira, Portugal: Evidence from vein spatial distributions, displacements and geometries[END_REF] proposed a model of vein opening involving episodic vein dilation induced by lithostatic fluid pressure conditions in compressive regime. To account for this feedback between permeability changes in the host rock domain and fluid overpressure conditions potentially induced by the expulsion of magmatic fluids, we evaluate the pore fluid factor at each time step:

(12) 𝜆 𝑣 = 𝑃 𝑓 𝜎 𝑣
where P f (Pa) is the fluid pressure, σ v (Pa) is the lithostatic pressure and λ v is the pore fluid factor, with λ v = 0.4 for hydrostatic fluid pressure conditions and λ v = 1 for lithostatic fluid pressure conditions. From values of λ v a factor of permeability change (f k ) was calculated assuming that it increases linearly with the pore fluid factor, with f k = 1 for λ v = 0.4 and f k = 100 for λ v = 1. The maximum value of f k was chosen according to the fluid pressure-dependent permeability relationships constrained from experimental studies [START_REF] Zhang | The influence of room temperature deformation on porosity and permeability in calcite aggregates[END_REF][START_REF] Faulkner | The effect of temperature, the nature of the pore fluid, and subyield differential stress on the permeability of phyllosilicate-rich fault gouge[END_REF][START_REF] Cox | Coupling between deformation, fluid pressures and fluid flow in ore-producing hydrothermal environments[END_REF][START_REF] Howald | Evidence for long-timescale (>103 years) changes in hydrothermal activity induced by seismic events[END_REF]. For values of λ v > 1 we assume that f k = 100. From this factor of permeability change, the permeability was calculated at each depth from its initial background value (Fig. 3b):

(13) 𝑘(𝑡) = 𝑘 0 .𝑓 𝑘
where k(t) (m²) is the permeability for a given depth and for a given fluid pressure conditions, k 0 (m²) is the background permeability at a given depth (Eq. 11) and f k is the factor of permeability change depending on the fluid pressure conditions.

Reaction-enhanced permeability during greisenization of granite

Porosity and permeability measurements performed by [START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF] on samples representative of different degrees of greisenization have emphasized a positive correlation between the porosity evolution and the alteration degree that follows an exponential law. To implement this experimental function of the porosity evolution in our numerical modeling, we have converted the alteration index (AI) in water/rock ratio (W/R) according to the results presented in the thermodynamic study of [START_REF] Reed | The Butte magmatic-hydrothermal system: One fluid yields all alteration and veins[END_REF]. This work uses batch equilibrium modeling to depict the stability of hydrothermal alteration assemblages encountered through the Butte granite during ore deposition. That means that W/R from this study is static, i.e., without fluid percolation through the porous matrix of the granite. However, as noted by [START_REF] Reed | Hydrothermal alteration and its relationship to ore fluid composition[END_REF], W/R can be used in an "approximate way" to describe the spatial and temporal evolution of alteration patterns in magmatichydrothermal systems: for example advanced argillic alteration is expected during "fresh" acidic fluids circulation and then large W/R for a "fluid-dominated" system, whereas a potassic-biotite alteration of granite is expected when fluids in lesser quantity have already reacted with rocks in a "rock-dominated" system with low W/R. In the Temperature vs. W/R diagram proposed by [START_REF] Reed | The Butte magmatic-hydrothermal system: One fluid yields all alteration and veins[END_REF] (Fig. 4a), the field of greisenization corresponds to the stability of muscovite (generated by acidic hydrolysis of feldspars and biotite from the granite) without recrystallizing feldspars. Therefore, from 550 to 300°C an optimum W/R value for greisenization reaction is close to 2. This does not imply that at W/R = 2 the granite is totally transformed into greisen and that feldspars and biotites no longer exist within the granitic matrix. It only means that greisenization is in progress and that only local fractions of feldspars and biotites are transformed into muscovite and quartz.

As above mentioned, the W/R in batch thermodynamic modeling is static. As porosity increases during greisenization [START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF] and that magmatic fluid supply is still effective up to the complete crystallization of the granitic magma (e.g., Lamy-Chappuis et al., 2020), we assume that "fresh" acidic fluids are renewed and recharge at each fluid flow pulses within the porous matrix. This corresponds to integrating into time successive fluid masses with respect to a constant rock mass to calculate a bulk time integrated W/R similar to that used in static thermodynamic batches (Fig. 4b). We have calculated the time-integrated water/rock ratios following the method proposed by Schardt and Large ( 2009):

(14) 𝑊/𝑅 = ∫ 𝑡 0 𝜌 𝑓 .𝑢.𝜑 𝜌 𝑟 .(1 -𝜑). 𝐿 .𝑑𝑡
where W/R is the time-integrated water/rock ratio, φ is the rock porosity, u is the fluid velocity, dt is the time difference between successive time step, ρ f is the fluid density, ρ r is the rock density and L the length scale of the defined mass rock, taken here as 1 m. Hence, W/R represents the total mass of fluid that has flowed through a defined mass of rock of scale L.

Therefore, from that optimum W/R value of 2 for greisenization reaction, we define a function describing the porosity evolution as a function of the W/R between 0 and 2 (Fig. 4c). This function conserves the same exponential evolution of the experimental law of [START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF]:

(15) 𝜑 (%) = 4.5 10 -1 .𝑒 (1.27. 𝑊/𝑅) At each time step of the calculation, the porosity of granite was calculated from this function for W/R ranging between 0 and 2 (Fig. 4b). For values of W/R ratio exceeding 2, we assume a fixed "greisen-like" porosity of 7%. Noteworthy, the W/R estimated to calculate the granite porosity is not dependent on the fluid source (magmatic from the intrusion vs. crustal-derived/meteoric within host rocks).

The permeability evolution related to the greisenization was assumed to vary with the porosity according to a cubic power law consistent with the experimental data described by [START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF] (Fig. 4d): ( 16)

𝑘 = 𝑘 0 ( 𝜑 𝜑 0 ) 3
where k 0 the initial granite permeability is 2.10 -20 m² and φ 0 the initial granite porosity is 0.4%.

Results

Following the approach described in section 3, we investigate through numerical modeling four different configurations which may be responsible for greisenization aiming at decipher specific effects of reactionenhanced permeability and magmatic fluid production. For these different scenarios, we first describe results obtained with static (constant) granite permeability. Then, permeability changes induced by greisenization are accounted for, following the approach described in section 3. The numerical modeling results obtained for ). For all models, during the cooling period we checked that phase separation (vapor) never occurred at any place (ESM5 and ESM6).

Models with static granite permeability

Without magmatic fluid production (M1)

This first model (Fig. 5) illustrates the development of a convective fluid flow pattern around the cooling granite. Despite convection clearly occurring, the small velocity values do not disturb significantly the isotherms, which remain parallel to the granite roof geometry. Temporal evolution shows that fluid loops, firstly focused just above the intrusion, progressively widen and propagate upward (black arrows in the left column of Fig. 5). This spatial propagation simply reflects the heat diffusion from granite, whose temperature equals 850°C at the initial time. The two convective cells appear to be centered above the cupola. This spatial control can be explained by strong contrasts induced by the granite geometry, such as large lateral temperature differences and high permeability contrast between the cupola and the host rock domain. The second column of Fig. 5 shows relatively low values of the time-integrated fluid flux. The maximum values are reached at time 50 kyr, and do not exceed 10 4.8 kg.m -2 . These maximum values appear localized along the granite-host rock contacts and above the cupola, which correspond to zones with the highest fluid velocities. In the third column, the fluid velocity curves show strong variations between the permeable host rock and the impermeable cupola, with the highest values (1.2 10 -9 m.s -1 ) at the incipient stage of the hydrothermal circulation.

With magmatic fluid production (M2)

The heat conduction and convection induced by the emplacement of granite intrusion lead to the progressive cooling of the magma chamber that triggers the production of magmatic fluids within the crystallized part of the magma chamber. This fluid production occurs approximately during the first 10 kyr and leads to the establishment of fluid overpressure conditions in the intrusive domain that reach lithostatic values (see white contours in Fig. 6, with λ v = 1). The pressure anomalies are localized along the granite roof during the first 5 kyr and extend progressively to the inner part of the granite in following the crystallization front through the magma chamber. The fluid overpressure conditions in the intrusive domain affect the fluid flow pattern, which is mainly marked by the expulsion of fluids from the upper part of the granitic intrusion toward the host rocks domain (Fig. 6). In this fluid flow pattern, the cupola appears to be a preferential zone of fluid expulsion. This expulsion of fluids in the host rock domain causes a moderate increase of the fluid pressure that enhances locally (i.e., mainly around the cupola) and moderately the permeability of the host rock domain (log of permeability change not exceeding 1.3). This slight permeability change suggests a limited effect of magmatic fluid production on the permeability evolution of the host rock domain.

Nevertheless, these zones of higher permeability are marked by fluid focusing and more elevated fluid velocities that locally reach 10 -8 m.s -1 . Moreover, extend of the area of permeability change in the host rock domain increase progressively over the first 8000 kyr. After 10 kyr, the granite is entirely crystallized, and the fluid production vanishes. However, lithostatic fluid pressure conditions remain within the intrusive body after 50 kyr (Fig. 6). This implies that the fluid produced during the granite crystallization is not efficiently extracted

Models considering reaction-enhanced permeability during the greisenization

Although the results described previously provide a useful reference framework for heat and fluid flow during the cooling stage of a granitic intrusion with a constant permeability, the effect of the metasomatic alteration on the granite permeability cannot be ignored. Here, the two following models include a timevarying permeability, induced by greisenization processes. In other words, equations ( 12) and ( 13) are implemented in the computation of the granite permeability for the two types of magmatic-hydrothermal systems presented previously.

Without magmatic fluid production (M3)

The left column of Figure 7 shows the same convective pattern as that observed in the case of constant granite permeability (Fig. 5), with small convective cells centered near the cupola and expanding with time.

The main difference consists in the progressive increase of the granite roof permeability. Indeed, as in the case of Figure 5, the fluid velocity is highest along the granite roof but in this case, permeability is increased (see grey zones in the left column) in the zones where the W/R ratio is high.

As a result, the greisenization front (red line in Figure 7) continuously deepens through the granite until the fluid velocity is not high enough to permit further fluid-rock interactions. Continuous development of this high permeability zone (the greisen) focuses fluid flow, as illustrated by a progressive increase of velocity (from 10 -9 m.s -1 at 5 kyr to 10 -8 at 50 kyr) across the cupola. This also causes elevated values of time-integrated fluid flux, up to 10 5.5 kg.m -² in the greisen, but also above the cupola over a large area (around 10 5 kg.m -²).

With magmatic fluid production (M4)

This last model (M4, Fig. 8) corresponds to the most elaborated one and includes both magmatic fluid production, permeability changes in the host rock domain due to potential fluid overpressure conditions, and permeability changes within the granite related to the greisenization processes.

The fluid flow pattern is firstly controlled by expulsion processes induced by magmatic fluid produced during the first 10 kyr (time step not shown in Figure 8). Then, due to complete granite crystallization, fluid production vanishes, and the convective process takes over the expulsion mechanism. As shown in the last time step, shallow convective cells are developed in the permeable upper part of the host rock domain.

Due to fluid production within the granite, the W/R ratio quickly increases in the upper part of the granite, leading to the formation of a permeable pathway along the granite roof (the greisen, shown in white in the left column). As in the previous case, the greisenization front deepens with time. It is noteworthy that the greisen thickness is slightly larger than in the previous case (M3).

In response to fluid expulsion, a fluid overpressured zone progressively develops above the granite roof, leading to a large zone of permeability increase that can reach two orders of magnitude. This zone quickly expands in the first 5 kyr, as shown by the vertical C-C' profile of permeability changes (right column). When the granite is completely crystallized (after ~ 10 kyr) fluid production ceases and consequently, the fluid pressure anomaly disappears.

The development of permeable zones in both the granite and the host rock is accompanied by high fluid velocity, reaching 10 -7 m.s -1 within and around the cupola (see right column). The isotherms, illustrated in the middle column (dashed white lines), are parallel to the pluton roof geometry in the first stages. In the last stage (50 kyr), the 200 to 400°C isotherms are stretched up by the hot fluids rising quickly above the cupola.

As far as time-integrated fluid flux is concerned, a similar description as in the M2 model can be made, but absolute values are one order of magnitude higher. A worthy difference, however, deals with the high value of the time-integrated fluid flux (up to 10 6 kg.m -²) within the upper part of the granite (greisen).

Discussion

This study deals with 2D numerical models investigating the influences of reaction-enhanced permeability related to greisenization processes, magmatic fluid production, and permeability changes in the host rock due to elevated fluid pressure conditions. Although some recent studies already explored the physical interactions between geological processes and permeability changes [START_REF] Weis | Porphyry-copper ore shells form at stable pressure-temperature fronts within dynamic fluid plumes[END_REF][START_REF] Mezri | New parametric implementation of metamorphic reactions limited by water content, impact on exhumation along detachment faults[END_REF][START_REF] Weis | The dynamic interplay between saline fluid flow and rock permeability in magmatic hydrothermal systems[END_REF][START_REF] Scott | Permeability Changes Resulting from Quartz Precipitation and Dissolution around Upper Crustal Intrusions[END_REF], our approach is also based on laboratory measurements of the evolution of porosity and permeability during greisenization. In the following, we compare numerical results from the four models M1 to M4 to discuss the effects of permeability changes on (i) expulsion of magmatic fluids, (ii) intensity of fluid fluxes and (iii) formation of massive greisens.

Effect of the cupola

In all cases (M1 to M4), the cupola strongly affects the fluid flow pattern by focusing both the convective upwelling and the expulsion of magmatic fluid. This can be explained by strong physical contrasts (permeability and temperature) between the cupola and its host rock. This result has already been described

by [START_REF] Eldursi | Numerical investigation of transient hydrothermal processes around intrusions: Heat-transfer and fluid-circulation controlled mineralization patterns[END_REF]2018). This fluid focusing by the cupola (from the cupola towards the host rock) has been independently deduced from a growth band study of tourmalines, which marks the incipient stage of the hydrothermal activity at Panasqueira [START_REF] Launay | Deciphering fluid flow at the magmatic-hydrothermal transition: A case study from the world-class Panasqueira W-Sn-(Cu) ore deposit (Portugal)[END_REF].

Besides focusing on fluid expulsion, plutonic apexes also localize thermo-mechanical instabilities (formation of fault and shear zones) [START_REF] Guillou-Frottier | The development and fracturing of plutonic apexes: implications for porphyry ore deposits[END_REF][START_REF] Gloaguen | Tracing the magmatic/hydrothermal transition in regional low-strain zones: The role of magma dynamics in strain localization at pluton roof, implications for intrusion-related gold deposits[END_REF]. The combination of these processes can favor the formation of hydrothermal systems, thus explaining the observed spatial relationship between cupolas and mineralized systems (e.g., [START_REF] Dilles | Metallogenesis of the Yerington batholith, Nevada[END_REF].

Interactions between dynamic permeabilities and extraction of magmatic fluids

The comparison between models M2 and M4 emphasizes that the permeability change induced by greisenization processes enhances and facilitates the expulsion of the magmatic fluids produced during granite crystallization. Figure 9 displays the time evolution of the pore fluid factor and the fluid flow velocity above and within the cupola. It shows that fluid overpressure conditions remain elevated in the granite (red curves)

-even after the end of fluid production -in model M2 (left column), while they drop quickly when the permeability change in granite is considered (M4). Moreover, fluid overpressure conditions above the granite (blue curves) are significantly more elevated in model M4 (λ v = 1) than in model M2 (λ v = 0.5). The highest fluid velocities both in the cupola and above the granite in model M4 are consistent with a more efficient fluid expulsion when the effects of greisenization on permeability are considered. This emphasizes that the increase of permeability induced by the greisenization causes the transfer of the overpressure generated by the production of magmatic fluids from the granite to the host rock domain above the granite roof. This efficient transfer of magmatic fluids across the greisen focalizes and stabilizes high fluid pressure conditions (λ v = 1) above the granite roof, which facilitates the development of a thicker band of high permeability in the host rock, as observed in Figure 8 at time 5 kyr. This zone of elevated permeability could correspond to the damage zone described by [START_REF] Foxford | Extensional failure and hydraulic valving at Minas da Panasqueira, Portugal: Evidence from vein spatial distributions, displacements and geometries[END_REF] (Fig. 1d), in which sub-horizontal veins could be preferentially opened in response to lithostatic fluid pressure conditions coupled with a compressive tectonic regime (i.e., horizontal effective major principal stress). Consequently, greisenization can constitute an important mechanism in the formation of large deposits by locating zones of higher permeability, which determines the locus of focused fluid expulsion, and thus enhances mass transfer within and above granitic intrusions. -16 to 10 -14.5 m²) in magmatic-hydrothermal systems where fluid fluxes reach 10 6 to10 7 kg/m² [START_REF] Gerdes | Convective fluid flow through heterogeneous country rocks during contact metamorphism[END_REF][START_REF] Cui | Controls of layered and transient permeability on fluid flow and thermal structure in contact metamorphic aureoles, with application to the Notch Peak aureole Utah[END_REF]. Our results suggest that dynamic permeability and greisenization processes could be important to promote and enhance the transport of a large amount of fluid and hence to form large deposits.

Influence of reaction-enhanced permeability on fluid flow and mass transfer

Implications for the formation of massive greisen

The interplays between greisenization and permeability have also a strong effect on the propagation of the greisenization front toward the deeper part of the intrusive body. Figure 11 shows that models M3 and M4 considering reaction-enhanced permeability during greisenization processes are characterized by a thicker greisen after 50 kyr of fluid flow than models (M1 and M2) with static granite permeability. Accordingly, the reaction-enhanced permeability seems to promote significantly the formation of massive greisen with thickness of about 140 m for the model M3 and 200 m for the model M4 (Fig. 11). Conversely, a static low permeability in granite allows only the development of a thin greisen with a thickness of about 20 m in the model M1 and 50 m in the model M2 (Fig. 11). Furthermore, model (M4) appears to be the most favorable for the development of massive greisen by vigorous fluid flow and self-propagating of pervasive alteration (greisenization). The positive feedback between greisenization, and permeability, and the magmatic fluid production could be important processes to form large massive greisen deposits (with disseminated mineralization) like Cinovec and East Kemptville, in which massive greisen can reach 200 to 300 m in thickness (Halter, 1996;[START_REF] Jarchovsky | The nature and genesis of greisen stocks at Krásno, Slavkovský les area western Bohemia, Czech Republic[END_REF][START_REF] Laznicka | Giant metallic deposits. Future sources of industrial metals[END_REF].

Limitations and further perspectives

Although results obtained in this study permit to discuss the potential role of interplays between reactionenhanced permeability and greisenization processes during the formation of vein and greisen deposits, we computed fluid flow on 2D sections with simplified geometries. Instead, 3D models would need to be considered in future works. Indeed, fluid flow patterns and fluxes described in this study should probably be different in 3D, but the interplays (physical processes) between dynamic permeability and greisenization processes and their consequences should not be affected by the geometry.

This study is being dedicated to feedback during reaction-enhanced permeabilities, the fracturing of solidified cupolas and granite roofs have not been considered as for example temperature-dependent permeabilities laws (considering effects of brittle-ductile transition on permeability) as explored by [START_REF] Hayba | Multiphase ground water flow near cooling plutons[END_REF]. More sophisticated magma chamber models accounting for hydraulic fracturing and capillary fracturing within mush have been published recently (e.g., [START_REF] Lamy-Chappuis | Mechanisms and patterns of magmatic fluid transport in cooling hydrous intrusions[END_REF]. They could be coupled with our approach and perhaps will lead to more vigorous systems with phase separation.

Furthermore, our models assume a single fluid phase (pure water). Recent studies [START_REF] Weis | Porphyry-copper ore shells form at stable pressure-temperature fronts within dynamic fluid plumes[END_REF][START_REF] Gruen | Hydrodynamic modeling of magmatic-hydrothermal activity at submarine arc volcanoes, with implications for ore formation[END_REF][START_REF] Weis | The dynamic interplay between saline fluid flow and rock permeability in magmatic hydrothermal systems[END_REF] have demonstrated that the salinity of fluid affects significantly (i) fluid flow patterns, (ii) fluid velocities, and especially (iii) fluid pressure conditions. In our models, the effects of greisenization and the opening of permeable zones above the granite should probably be enhanced in the case of a multiphasic fluid. Indeed, the potential increase in fluid volume induced by the separation of a vapor phase can increase drastically fluid pressure conditions favoring the opening of permeable zones [START_REF] Weis | The dynamic interplay between saline fluid flow and rock permeability in magmatic hydrothermal systems[END_REF].

The stress state is another parameter not considered in this study and that could also enhance and promote the opening of permeable zones [START_REF] Zhang | The influence of room temperature deformation on porosity and permeability in calcite aggregates[END_REF][START_REF] Cox | The application of failure mode diagrams for exploring the roles of fluid pressure and stress states in controlling styles of fracture-controlled permeability enhancement in faults and shear zones[END_REF][START_REF] Weis | The dynamic interplay between saline fluid flow and rock permeability in magmatic hydrothermal systems[END_REF].

Results of modeling described in this study deal with bulk pervasive fluid flow around and within granite intrusion, i.e., at the thermal aureole scale. At this scale, veins are not identifiable, and localized fluid flows around and within individual opening fractures could not be addressed in this study.

Besides, metamorphic reactions related to the thermal metamorphism induced by the granite emplacement can lead to the production of metamorphic fluids (dehydration reactions) that could constitute an additional source of fluids during the incipient stage of the magmatic-hydrothermal systems [START_REF] Connolly | Fluid and enthalpy production during regional metamorphism[END_REF]Hanson, 1992;[START_REF] Cui | Controls of layered and transient permeability on fluid flow and thermal structure in contact metamorphic aureoles, with application to the Notch Peak aureole Utah[END_REF]. This additional source of fluids can promote fluid overpressure conditions and enhance the opening of permeable zones above the granite intrusion. Further tests need to be performed to constrain the influence of thermal metamorphism on fluid flow dynamics. A recent study [START_REF] Scott | Permeability Changes Resulting from Quartz Precipitation and Dissolution around Upper Crustal Intrusions[END_REF] has demonstrated that quartz precipitation and dissolution can lead to significant permeability changes. In our study, we have not considered the potential precipitation of quartz in the neoformed porosity of greisen that could partially seal pathways and inhibit the feedback between permeability and greisenization processes.

Additionally, we simplified the magmatic-hydrothermal system by a single stage of pluton emplacement and its subsequent cooling. The concept of incremental emplacement of granitic bodies [START_REF] Dilles | The petrology of the Yerington Batholith, Nevada: evidence for the evolution of porphyry copper ore fluids[END_REF][START_REF] Sillitoe | Porphyry Copper Systems[END_REF][START_REF] Blundy | Crustal magmatic systems from the perspective of heat transfer[END_REF] is thus not considered, but our results show that geological processes investigated in this study occur in a very short time span. Consequently, our hypothesis can be considered as representative of one pulse of a long time magmatic-hydrothermal history.

The emplacement depth of the intrusion is a key-controlling factor in the hydrodynamics and thermodynamics of such intrusion-related mineralized systems (e.g., [START_REF] Eldursi | Numerical investigation of transient hydrothermal processes around intrusions: Heat-transfer and fluid-circulation controlled mineralization patterns[END_REF][START_REF] Chi | Hydrodynamic links between shallow and deep mineralization systems and implications for deep mineral exploration[END_REF].

Numerical investigations accounting for emplacement depth coupled with petro-physical feedback we present here are perspectives to be explored in the close future. Despite these numerous limitations and future perspectives, the present study provides an additional contribution describing the interplays between hydrothermal alteration, dynamic permeability, the expulsion of fluid, and formation of overpressured permeable zones in magmatic-hydrothermal systems. These complex interactions can constitute important mechanisms to enhance and promote mass transfer to form large vein and greisen ore mineralized systems.

Finally, this contribution demonstrates that feedbacks between fluid-rock reactions, porosity, and permeability changes affect significantly fluid flow and need to be considered in future numerical modeling studies to provide more accurate predictive models.

Table 1 Initial physical properties of rocks and fluid used for models. Eq(n) refers to equations presented in the text.

Figure 1 (a) Tectonic map displaying the tectonostratigraphic domains composing the Iberian massif (modified from [START_REF] Julivert | Mapa Tectonico de la Peninsula Iberica y Baleares E[END_REF]. CZ: Cantabrian Zone, WALZ: West Asturian-Leonese Zone, GTOMZ: Galicia Tràs-os-Montes Zone, CIZ: Central Iberian Zone, OMZ: Ossa-Morena Zone and SPZ: South Portuguese Zone. (b) Regional geological map of the W-Sn-(Cu) Panasqueira ore deposit, in which the lithogeochemical anomalies in W are represented (compiled from data of [START_REF] Oosterom | Lithogeochemical studies of aureoles around the panasqueira tin-tungsten deposit[END_REF]. (c) The (A-A') geological cross section displays the spatial relationships between the mineralized vein system, the greisen and the Panasqueira granite. (d) Cumulative vein thickness along the vertical sections of drill holes showing the vertical distribution of the mineralized vein network of Panasqueira (n refers to the number of veins crosscut by the drill holes) (drill holes data compiled from [START_REF] Foxford | Extensional failure and hydraulic valving at Minas da Panasqueira, Portugal: Evidence from vein spatial distributions, displacements and geometries[END_REF][START_REF] Wheeler | Technical report on the mineral resources and reserves of the Panasqueira mine, Portugal[END_REF]. [START_REF] Manning | Permeability of the continental crust: implications of geothermal data and metamorphic systems[END_REF]. Permeability increase induced by fluid overpressure conditions (λ v > 0.4) are also displayed. For a given depth, the permeability increases by 2 orders of magnitude from the background value when the fluid pressure conditions are lithostatic (λ v = 1). (b) Permeability evolution at 2 km of depth as a function of the pore fluid factor (λ v ), shown here with a logarithmic scale. We assume that the permeability increases linearly with the pore fluid factor until lithostatic pressure conditions (λ v = 1). ). We assume that the W/R ratio of 2 is an optimum value to achieve efficient greisenization of the granite (see explanation in section 3.6.3). The red line indicates the emplacement of the profile along which, the time-evolution of the greisen thickness has been determined.

  these different configurations are displayed in Figures 5 to 8 as snapshots of different times illustrating the time evolution of the fluid flow patterns and the key processes involved during fluid flow. For each model, one figure presents three-time steps of fluid flow patterns and permeability changes (left column), thermal field and time-integrated fluid flux (middle), and cross sections of some hydrodynamics parameters through the cupola (horizontal B-B' profile) and along the C-C' vertical profile (right column). Additional figures displaying for the same three time steps the pressure and temperature conditions, as well as the fluid density and the porosity of granite for each model are available in electronic supplementary materials (ESM1 to ESM4

  from the intrusive body due to the low permeability of the granite domain. This fluid pressure anomaly affects strongly the convective fluid flow processes that dominate the fluid flow pattern at 50 kyr. Indeed, convective cells are marked by an important upward fluid flow above the cupola that is probably induced by the slow and continuous expulsion of fluids from the granite domain.During the first 10 kyr, the highest values of time integrated fluid flux are mainly localized above the granite roof and the cupola due to the expulsion of fluids produced during the granite crystallization. After 50 kyr of fluid flow, the time integrated fluid flux reaches values of 10 5 kg.m -2 . These maximum values are mainly concentrated above the cupola that corresponds to zones of fluid focalization and highest fluid velocities in response to the permeability increase observed during the first 10 kyr. The presence of a ~200 m thick zone with relatively high time-integrated fluid flux (light blue zone at 50 kyr, 10 3.75 kg.m -2 ) in the upper part of the granite is consistent with the preferential release of magmatic fluid from the granite roof.

Figure 10

 10 Figure10shows the time evolution of fluid fluxes for the four models, computed in a small area above the cupola (red area). Models M1 and M3 (without fluid production) are shown in grey. As expected, the model with reaction-enhanced permeability in granite (M3) is marked by a higher fluid flux, with a maximum value of ~10 -2 kg.s -1 reached at 50 kyr. Maximum fluid fluxes are achieved for M2 and M4 (black curves, with fluid production) at the end of the fluid production stage (~10 kyr). After this stage, a regular decrease in the fluid flux illustrates the progressive cooling of the magmatic-hydrothermal system. The account of (i) a dynamic permeability both in the granite and in the host rocks and (ii) a magmatic fluid source, leads to the highest fluid flux (3.10 -1 kg.s -1 , equivalent to ~ 10 Mt.kyr -1 ) (model M4, dashed black line in Figure11a). Due to the development of high permeability zones in the upper part of granitic intrusion, the values of timeintegrated fluid flux in models (M3 and M4) considering the reaction-enhanced permeability in granite (Figs.7 and 8) are higher by approximately one order of magnitude than values obtained in their equivalent models considering static granite permeability (Figs.5 and 6). It appears that model M4 is the most favorable for transferring a large amount of fluid. As recently confirmed by[START_REF] Cernuschi | Rapid formation of porphyry copper deposits evidenced by diffusion of oxygen and titanium in quartz[END_REF], a large amount of focused mineralizing fluids is required to form large magmatic-hydrothermal ore deposits. Furthermore, the highest values of time integrated fluid flux obtained in the model M4 (with reaction-enhanced permeability and fluid production) are similar to those described in numerical studies involving high values of static permeability (10 -16 to 10 -14.5 m²) in magmatic-hydrothermal systems where fluid fluxes reach 10 6 to10 7 kg/m²

Figure 2

 2 Figure2Meshing and model geometry with a laccolith magma chamber of about 6 km of width for 2 km of height and a cupola at 2.7 km of depth. This geometry corresponds to the inferred geometry of the Panasqueira granite. The boundary and the initial conditions used for the modeling are also displayed. Cross-sections B-B' and C-C' refer to horizontal and vertical profiles used in Figures 6 to 9.

Figure 3

 3 Figure 3 Model of permeability in the host rock domain describing the permeability increase induced by fluid overpressure conditions. (a) Permeability profile as a function of depth (adapted from Ingebristen and[START_REF] Manning | Permeability of the continental crust: implications of geothermal data and metamorphic systems[END_REF]. Permeability increase induced by fluid overpressure conditions (λ v > 0.4) are also displayed. For a given depth, the permeability increases by 2 orders of magnitude from the background value when the fluid pressure conditions are lithostatic (λ v = 1). (b) Permeability evolution at 2 km of depth as a function of the pore fluid factor (λ v ), shown here with a logarithmic scale. We assume that the permeability increases linearly with the pore fluid factor until lithostatic pressure conditions (λ v = 1).

Figure 4

 4 Figure4(a) Alteration mineral phase diagram displaying the stability fields of minerals for the Butte granite as a function of the water/rock ratio and temperature conditions (modified from[START_REF] Reed | The Butte magmatic-hydrothermal system: One fluid yields all alteration and veins[END_REF]. The red dashed line marks the total breakdown of primary magmatic minerals (albite, k-feldspars and biotite) and the filled area corresponds to the conditions favourable for greisenization. Mineral abbreviations: Ab = Albite, Bt= Biotite, Kfs= k-feldspar, Ms= Muscovite. (b) Conceptual model showing coupling between the W/R ratio and reaction-enhanced permeability in the granite during greisenization. (c) Model of porosity evolution in greisen as a function of the water/rock (W/R) ratio and the alteration index (AI = 100.(Al 2 O 3 /(Al 2 O 3 +K 2 O+Na 2 O))). The conversion of the alteration index (AI) in water rock (W/R) ratio is detailed in the text. (d) Model of permeability evolution in greisen as a function of the porosity generated during the greisenization. The red curves correspond to functions (Eq. (15) and Eq. (16)) used in numerical modeling to describe the porosity and the permeability evolution during the greisenization (Experimental data after[START_REF] Launay | Dynamic permeability related to greisenization in Sn-W ore deposits: Quantitative petrophysical and experimental evidence[END_REF]).

Figure 5

 5 Figure 5 Results of the time evolution of the fluid flow patterns and the time-integrated fluid flux obtained for modeling with static granite permeability and without magmatic fluid production (M1). After 50 kyrs a near steady state conditions are achieved.

Figure 6

 6 Figure 6Results of the time evolution of the fluid flow patterns and the time-integrated fluid flux obtained for modeling with static granite permeability considering magmatic fluid produced during the granite crystallization (M2). Fluid flow is mainly controlled by expulsion processes during the first 10 kyrs and by convective circulation after 50 kyrs, (the localization of the vertical and the horizontal profiles are given in Figure2).

Figure 7

 7 Figure 7 Results of the time evolution of the fluid flow patterns and the time-integrated fluid flux obtained for modeling considering reaction-enhanced permeability during the greisenization process in granite and without fluid production (M3). Fluid flow is mainly controlled by convective processes.

Figure 8

 8 Figure 8 Results of the time evolution of the fluid flow patterns and the time-integrated fluid flux obtained for modeling considering both reaction-enhanced permeability and magmatic fluid production during the

Figure 9

 9 Figure 9 Time-evolution of pore fluid factor and fluid flow velocity above and within the cupola (emplacements indicated by the blue and the red stars) in cases of magmatic fluid production without (left) and with (right) reaction-enhanced permeability in granite. Light grey areas correspond to the timespan of fluid production and dark grey areas highlight the sharp increase of fluid velocity and overpressure in the M4 model.

Figure 10

 10 Figure 10 Effects of the reaction-enhanced permeability and the magmatic fluid production on the time evolution of fluid flux above the cupola (red area).

Figure 11

 11 Figure 11Time-evolution of the greisen thickness resulting from fluid-rock interactions during fluid flow (for M1 to M4). We assume that the W/R ratio of 2 is an optimum value to achieve efficient greisenization of the granite (see explanation in section 3.6.3). The red line indicates the emplacement of the profile along which, the time-evolution of the greisen thickness has been determined.

Acknowledgments

This research took part in the ERAMIN project "New Ores" (ANR-14-EMIN-0001) and was financially supported by the French Geological Survey (BRGM), the Region Centre, and the Labex Voltaire (ANR-10-LABX-100-01). We are grateful to anonymous reviewers for their constructive remarks, which have greatly improved this manuscript.

Authorship contributions statement:

Gaëtan Launay: Conceptualization, methodology, numerical modeling and formal analysis of results, visualization and construction of figures, writing of the original draft, and editing. Eq. ( 10) kg of fluid/m 3 of rock/s (kg/m 3 /s)

Highlights

-Reaction-enhanced permeability during greisenization promotes the expulsion of magmatic fluids.

-The efficient expulsion of magmatic fluids during greisenization stabilizes fluid overpressures responsible for the opening of vein swarm within the host rocks.

-Interplays between permeability and greisenization promote and enhance the transport of a large amount of fluid.

-Reaction-enhanced permeability is a key mechanism for formation of large vein and greisen Sn-W deposits.