Regional thermo-rheological field related to granite emplacement in the upper crust: implications for the Larderello area (Tuscany, Italy)

ABSTRACT

We modelled thermo-rheological perturbations, related to the emplacement of a magmatic body in the upper crust. This approach was considered relevant for the areas characterized by elevated surface heat flow and chiefly for the geothermal fields. The numerical conductive thermal model applied to the Larderello geothermal area in Tuscany, allowed to constrain size, depth and timing of emplacement of the pluton. We inferred that the emplacement of a magmatic body, at a minimum depth of 3 km, having a horizontal extension of 14 km and a maximum thickness of 8 km, can reasonably reproduce the observed regional surface heat flow anomaly of the Larderello area, when 300 (± 100) kyr are elapsed from the magma emplacement. Even assuming an incremental growth, the first magma injection should not be older than 1 ± 0.3 Ma.

Results of the thermal model were used to set up a rheological model and to simulate the drifting of the brittle-ductile transition during the cooling of the pluton. A comparison with the K-horizon profile, a prominent seismic reflector in the Larderello area, was then performed. It was found that the K-horizon approximately corresponds with the pluton roof and with the current location of the brittle-ductile transition.

KEYWORDS:

• Numerical thermal model
• heat flow anomalies
• brittle-ductile transition (BD¹)
• K-horizon

1. Introduction

The surface heat flow is a key parameter for understanding the thermal conditions and processes acting at depth (Ranalli & Rybach, 2005). It results from various contributions such as (i) the heat flowing from the mantle into the crust, (ii) the radiogenic heat produced within the continental crust, (iii) tectonothermal events such as lithospheric thinning or thickening, (iv) surface and shallow crustal levels processes such as sedimentation, erosion, uplift, subsidence, paleoclimate T changes, fluid circulation and, finally, (v) the emplacement of igneous bodies in the upper crust (e.g., Della Vedova, Bellani, Pellis, & Squarci, 2001). The relationship between high surface heat flow areas and magma emplacement in the continental crust is well documented in the literature (Allis et al., 1995; Baldi et al., 1994; Cathles, Erendi, & Barrie, 1997; Gianelli, Manzella, & Puxeddu, 1997; Gupta, 1981; Reiter, Chamberlain, & Love, 2010; Stimac, Goff, & Wohletz, 1997, 2001). The increase of surface heat flow is particularly elevated when magma emplacement takes place at shallow depth (e.g., Stimac, Goff, & Goff, 2015). To quantify this effect, thermal modelling has been performed by considering both (i) conductive (Allis et al., 1995; Dalrymple et al., 1999; Rikitake, 1995; Stimac et al., 2001) and (ii) convective heat transfer (e.g., Dalrymple et al., 1999; Stimac et al., 2001). Norton & Knight (1977) simulated cooling of a magmatic body emplaced at shallow depth and demonstrated that the resulting temperature field obtained with the two approaches is significantly different in zones located above the pluton roof but is comparable in zones located at the margin or in the surroundings of the pluton. On this basis, it could be reasonable to apply conductive models to geothermal areas if the objective is circumscribed to the simulation of the long wavelength thermal anomaly related to magma emplacement (Dalrymple et al., 1999). By reproducing the regional thermal anomaly is then possible to constrain size, depth and age of emplacement of the magmatic body. In this paper, we applied this approach to the well-known high heat flow area of Larderello in Tuscany. Here, the magmatic source, invoked in all studies, is still poorly defined, despite the importance of the related geothermal field, the most ancient in the world (e.g. Saccorotti et al., 2014).

A debatable feature of the Larderello area is the prominent seismic reflector, called K-horizon (e.g. Cameli, Dini, & Liotta, 1993), whose origin has not understood yet. For this reason, once defined the time-dependent thermal profiles of the area, a rheological model has been developed with the aim to provide a contribution to the interpretation of the K-horizon.

2. Relevant features of the Larderello geothermal area

The Larderello geothermal area (Figure 1) is located in the inner part of the Northern Apennines chain (southern Tuscany). This area, since early-middle Miocene, has been affected by extensional tectonics that produced a relevant lithospheric thinning and favoured the development of Neogene basins with dominant NW-SE orientation. Current thickness of crust and lithosphere are estimated at c. 23 km (e.g., Ponziani et al., 1995) and c. 40 km (e.g., Calcagnile & Panza, 1980), respectively. Moreover, this area was affected by magmatism, that migrated progressively eastwards, from Corsica (Sisco, 14 Ma) to the inner Tuscany (Mt. Amiata, 0.2 Ma), being responsible for the genesis of the acid intrusive bodies of Montecristo, Giglio and Elba Island and of the volcanic products of Capraia and Mt. Amiata (e.g., Dini, Gianelli, Puxeddu, & Ruggieri, 2005). According to Peccerillo (2005), the origin of silicic magmatism in Tuscany is attributed to crustal anatexis.

Figure 1. Regional heat flow contour lines in southern Tuscany. Equidistance is 50 mWm⁻² (modified after Batini, Brogi, Lazzarotto, Liotta, & Pandeli, 2003).

In the Larderello area, the local heat flow anomaly, dominated by convection, reaches a peak value in the order of 1000 mWm⁻² (e.g., Bellani, Brogi, Lazzarotto, Liotta, & Ranalli, 2004). Instead, the regional long wavelength heat flow anomaly (Figure 1), dominated by conduction, has a maximum in the order of 150–200 mWm⁻² (e.g., Bellani, Magro, Brogi, Lazzarotto, & Liotta, 2005), well above the average value for a continental area of 65 mWm⁻² (Pollack, Hurter, & Johnson, 1993). Hydrothermal circulation is responsible for the short wavelength thermal anomaly (Baldi et al., 1994; Magro, Bellani, & Della Vedova, 2009), whereas extensional tectonics and emplacement of magmas in the upper crust have the main influence on the long wavelength thermal anomaly (Baldi et al., 1994; Bellani et al., 2004). Both convective (Della Vedova, Vecellio, Bellani, & Tinivella 2008; Ebigbo et al., 2016; Magro et al., 2009; Romagnoli, Arias, Barelli, Cei, & Casini, 2010) and conductive (Calore, Celati, Gianelli, Norton, & Squarci, 1981; Mongelli, Palumbo, Puxeddu, Villa, & Zito, 1998) models have been proposed for this area.

This paper is focused on the interpretation of the regional heat flow anomaly and its relationship with the magmatic heat source. In the Larderello geothermal field, the heat source is attributed to shallow igneous intrusions, belonging to the Tuscan Magmatic Province (Gianelli, 2008, and references therein). The presence of a magmatic body is confirmed by most of the deep wells that reached granite intrusions (Dini et al., 2005; Villa & Puxeddu, 1994). However, despite the numberless studies on the Larderello area, there are still different opinions about the size, depth and age of emplacement of the magmatic heat sources underneath the geothermal field. Del Moro, Puxeddu, Radicati Di Brozolo, & Villa (1982), attributed the origin of the geothermal field to the intrusion of a granitic body having an age older than 3.7 Ma and a depth of emplacement lower than 9 km. According to Villa & Puxeddu (1994), the emplacement of the magma took place at depths ranging from 4 to 9 km and occurred at 4 Ma. Afterwards, the granitic body underwent a monotonic cooling, without evidence of ascent of multiple magmatic pulses. The still active geothermal anomalies are attributed to the predominance of slow conduction over fast convection heat transfer (Villa & Puxeddu, 1994). Mongelli et al. (1998) built up a conductive model based on the cooling of a huge magmatic body in order to reproduce the temperatures observed in the wells of the geothermal field. Their results suggest that the last magma injection occurred at 0.3 Myr. According to Dini et al. (2005), multiple intrusions of granitic magmas emplaced between 3.8 and 1.3 Myr, at depths ranging from 3 to 5 km. The geological cross-sections of Bertini, Casini, Gianelli, & Pandeli (2006) show a convex-upward shape of a Quaternary granite, up to 20 km wide, emplaced at 3–4 km depth. Rossetti et al. (2008), suggested that the pluton, in the Larderello area, is the result of the coalescence of multiple intrusions. During the growth, roof lifting provided the room for magma emplacement. More recently, Spinelli et al. (2015) inferred that dated intrusions (see Table 1) cannot be responsible for the active geothermal system, considering the short time (in the order of 10–100 kyr) sufficient to reach thermal equilibrium between the magmatic bodies and wall rocks. Santilano et al. (2015) sustained that the shallow magmatic intrusion, underneath the Larderello geothermal system, should be more recent than 1 Myr. The same authors, in accordance with Spinelli et al. (2015), attributed the thermal anomaly to the shallow emplacement of a very recent granitic magma body.

Table 1. Depth and age of the magmatic rocks cored in the Larderello area.

Sample	Depth (m)	Age (Ma)	Method	Reference
Monteverdi 7	3483	3.8 ± 0.1	K/Ar (biotite)	Villa, Gianelli, Puxeddu, Bertini, & Pandeli (1987)
Montieri	2200	3.1 ± 0.1	K/Ar (biotite)	Villa, Ruggieri, Puxeddu, & Bertini (2006)
Montieri	2670	3.7 ± 0.1	K/Ar (biotite)	Villa et al. (2006)
Montieri	3400	2.9 ± 0.1	K/Ar (biotite)	Villa et al. (2006)
Sasso 22	3529	3.06 ± 0.07	K/Ar (biotite)	Villa & Puxeddu (1994)
Sasso 22	4028	2.88 ± 0.05	K/Ar (biotite)	Villa & Puxeddu (1994)
VC-11	2946	2.89 ± 0.05	K/Ar (biotite)	Batini et al. (1985a)
Radicondoli 26B	4120	2.47 ± 0.02	^40Ar/^39Ar (muscovite)	Villa, Ruggieri, & Puxeddu (1997)
Radicondoli 29	4300, 4548, 4822	2.25 ± 0.04	^40Ar/^39Ar (biotite)	Gianelli & Laurenzi (2001)
San Pompeo 2	2963	1.61 ± 0.12	K/Ar (biotite)	Batini et al. (1985a)
Carboli C bis	4200	1.23 ± 0.13	^40Ar/^39Ar (muscovite)	Villa et al. (2001)
Carboli C bis	4304	1.27 ± 0.22	^40Ar/^39Ar (muscovite)	Villa et al. (2001)

Regarding geophysical studies, seismic data revealed the presence of three low velocity zones below the Larderello area: two anomalies at depths of 4 and 7 km, with a width ranging from 5 to 10 km, and a third anomaly located at depths between 7 and 20 km (Gianelli et al., 1997, and references therein). Also, magnetotelluric surveys indicated the presence of a large conductive body at a depth of 5–6 km (Fiordelisi et al., 1995). Mongelli et al. (1998), on the basis of tomographic studies of Foley, Toksoz, & Batini (1992), suggested a cylindrical shape for the magmatic body, a width of 15–20 km and a height of 40 km.

Summing up, from the inspection of the literature on the Larderello area, it emerges that the magmatic heat source responsible for the thermal anomaly has a size not well defined and a depth of emplacement ranging between 3 and 9 km. Regarding the age, it comes up a dichotomy between geochronological dating (4–1.3 Myr) and thermal model results (< 1–0.3 Myr). In addition, both a single and a multiple intrusion of granitic magma have been proposed.

Another interesting feature of this area regards the so called ‘K-horizon’, a prominent seismic reflector (Figure 2), regionally located at depths of 8–10 km, and showing an evident culmination in the Larderello and Mt. Amiata geothermal areas, where it reaches depths of 3–6 km (e.g. Cameli, Dini, & Liotta, 1998) (Figure 3). The interpretation of this reflector is still controversial. For some authors, it represents the roof of a kinematically active shear zone, located at the top of the BD¹ transition (Brogi, Lazzarotto, Liotta, Nicolich, & Ranalli, 2003; Cameli et al., 1998; Liotta & Ranalli, 1999). The strong contrast in acoustic impedance, causing the high amplitudes, could be explained by the presence of shear zones with mylonitic rocks (Batini, Duprat, & Nicolich, 1985b; Cameli et al., 1993) and very high pore fluid pressure (Brogi et al., 2003; Liotta & Ranalli, 1999). Also, based on deep boreholes, the depth of the K-horizon can be correlated with an isotherm in the range of 450 ± 50°C (Cameli et al., 1998; Liotta & Ranalli, 1999). According to other authors, the K-horizon could be related to the presence of a granitic body. In particular, Gianelli et al. (1988), suggested that this reflector, in the Mt. Amiata geothermal area, could correspond to a fractured zone, with thermo-metamorphic minerals and fluids, at the top of an igneous body. According to Bertini et al. (2006), the K-horizon, in the Larderello area, corresponds to the roof of a Quaternary granitic intrusion and related fractured contact rocks permeated by a supercritical fluid.

Figure 2. Line-drawing showing the profile of the K-horizon as deduced from the seismic section in the Larderello-Travale geothermal area (as reported in Batini et al., 2003).

Figure 3. Contour lines in kilometers (equidistance: 0.5 km) of the K-horizon depth (modified after Batini et al., 2003). The blue line is the trace of the K-horizon profile shown in Fig. 14–15.

Figure 14. Rheological cross-sections at 300 kyr after the magma emplacement. The colour scale denotes Log₁₀[strength, MPa]. The black line shows the location of the BD¹ transition. The blue dots indicate the depth of the K-horizon along the profile, represented by the blue line in Figure 3.

Figure 15. Temperature (°C) cross-section at 300 kyr after the magma emplacement. The white dashed line represents the boundary of the magmatic body. The black line shows the location of the BD¹ transition obtained from the rheological calculations. The blue dots indicate the depth of the K-horizon along the profile represented by the blue line in Figure 3.

3. Thermal model

A method to obtain information on the thermal and geometrical properties of the heat source in the Larderello area is the reproduction of the observed regional thermal anomaly and the comparison between model estimates and observed data (e.g. de Lorenzo, Di Renzo, Civetta, D’Antonio, & Gasparini, 2006). In this way, we can gauge extent, depth and age of magmatic body emplacement. We have constrained these parameters by using a trial and error approach, in order to fit the observed regional heat flow anomaly (150–200 mWm⁻²) along a NW-SE profile (Figure 4). For this purpose, we set up a 2D numerical model to solve the following heat equation:

(1)

Figure 4. Heat flow contour map (mWm⁻²) in the Larderello geothermal field, modified after Magro et al. (2009). The red line is the trace of the surface heat flow profile shown in Figure 8 (red dots).

Figure 8. Comparison between observed (red dots) and model (black line) values at 300 kyr of surface heat flow along the profile indicated by the red line in Figure 4. Heat flow range in the 200–400 kyr time interval is indicated in grey.

where is temperature; is time; and are horizontal and vertical coordinates; is thermal conductivity; is rock density, is specific heat at constant pressure and is radiogenic heat production. The numerical solution to equation (1) has been obtained adopting the explicit solution to the finite difference problem (e.g., Anderson, 1995), using first order forward difference for the computation of the time derivative and second order central differences for the computation of the spatial derivatives. This approach allows the computation of the temperature at each grid point (i∆x, j∆z) (i = 1,…N; j = 1,…M), at the time step (k + 1) ∆t, if the temperature in the adjacent points, at the time k∆t is known. The time marching scheme is the following:

(2)

where is the thermal diffusivity, and are, respectively, the temperature and the internal heat generation in the grid points () at time step . For the explicit approach, the value of is chosen according to the following stability criterion (e.g., Becker & Kaus, 2016):

(3)

where and indicate the node spacing in spatial direction and , respectively. In this case, has been assumed.

The thermal conductivity is assumed equal to 2.0 W m⁻¹°C⁻¹. This represents the mean value of thermal conductivity in the earth’s crust (e.g., Kappelmeyer & Hänel, 1974). Even, considering a thermal dependence of thermal conductivity, the results in terms of size, depth and age of magmatic body are substantially unchanged. For this reason, we adopted the simplest configuration, assuming a constant thermal conductivity.

The initial geotherm has been computed by assuming T = 10°C at the surface (z = 0 km) and T = 900°C at the bottom of the crust (z = 23 km). The latter value is required for significant biotite dehydration melting in the lower crust (e.g., Bucher & Grapes, 2011). It was also assumed that the radiogenic heat production decreases exponentially with depth, according to the following expression (e.g., Lachenbruch, 1970):

(4)

where is the surface (z = 0) radioactive heat generation (3 μWm⁻³) and is the scale length of radioactivity (15 km) (e.g., Turcotte & Schubert, 2014). In this way, the initial surface heat flow is equal to 100 mWm⁻² (Figure 1).

The contribution to the temperature field due to the latent heat of crystallization has been computed by using the approach proposed by Spear (1993), according to the following equation:

(5)

where is the effective specific heat (2500 J kg⁻¹°C⁻¹), is the true specific heat (1000 Jkg⁻¹°C⁻¹), is the latent heat of crystallization equal to 300 kJ kg⁻¹ and is the temperature interval of magma crystallization (850–650°C) (e.g. Spear, 1993). A typical density value of 2750 kg m⁻³ has been assumed for the rocks in the crust.

A simplified hexagonal shape has been adopted in order to approximate the geometry of a lenticular magmatic body (Figure 5). This shape may be in agreement with a geological cross-section of Bertini et al. (2006), at least for the roof of the intrusion. Since we have not insights about the lower part of the magmatic body, a symmetrical shape has been considered, that may be produced by mechanisms of roof lifting and floor depression (Cruden, 1998).

Figure 5. Schematic figure of the geometric shape used to approximate a lenticular magmatic body. The black arrowed lines indicate the size (w = width, h₁ = lateral thickness, h₂ = central thickness) and the depth (d) of the pluton constrained, from the surface heat flow data, through a trial and error approach (see later).

In order to model this shape, the discretization step is computed, after setting discretization step and hexagon dimension, following this expression:

(6)

where is the slope of one side of the hexagon. In this way, the borders of the body cross the grid points, avoiding the calculation of their position through interpolation. In this case, = 0.5 km and = 0.14 km.

4. Rheological model

In order to obtain 2D strength profiles of the upper crust (down to a depth of 15 km) and to define the boundary between the brittle and ductile domains, we set up a rheological model following principles and methods described by Ranalli (1995). In addition, in order to provide a contribution to the interpretation of the K-horizon, we compared the depth of this seismic reflector with depth of BD¹ transition.

The brittle strength has been computed by using the Coulomb–Navier shear failure criterion (e.g. Liotta & Ranalli, 1999):

(7)

where and are maximum and minimum principal stresses, is the average rock density above depth z, is the gravity, is the pore fluid factor (ratio of fluid pressure to lithostatic pressure) and is a dimensionless parameter depending on the predominant tectonic regime and on the friction coefficient of rocks. We assumed an extensional regime (= 0.75; see Ranalli, 1995), equal to 2750 kg m⁻³ and

0.85 << 0.99 (see Liotta & Ranalli, 1999). In particular, in the country rocks located below and aside the magmatic body, increases linearly from 0.85 (at z = 0 km) to 0.99 according to the following expression:

(8)

thus, is equal to 0.99 at z = 5 km and this value is kept constant down to 15 km. Instead, for the country rocks located above the magmatic body, we adopted an expression accounting for the dependence of on both z and x. Therefore, increases from 0.85 at the surface to 0.99 at the interface between the host rocks and the magmatic body, according to this equation:

(9)

This takes into account the increase with the approach to the pluton. Finally, a equal to 0.99 has been adopted for the magmatic body (Figure 6).

Figure 6. Schematic figure of the pore fluid factor () variation above the roof of the pluton (dashed white line).

The ductile strength has been computed using the power-law creep equation (e.g. Ranalli, 1995):

(10)

Where is the strain rate, and are the maximum and minimum principal stresses,, and are the creep parameters, is the gas constant and is the absolute temperature. This equation can be expressed in terms of stress as:

(11)

The 2D temperature field obtained from the previous thermal model is used to estimate the ductile strength in the upper crust. We adopted a schematic lithological section (Figure 7), on the basis of the geological cross-sections through the Larderello geothermal area (e.g. Brogi & Liotta, 2006), assuming the flow law parameters listed in Table 2. A critical point concerns the choice of the strain rate parameter. In fact, considerably different values have been assumed or estimated in the past (e.g., Acocella & Rossetti, 2002; Dallmeyer & Liotta, 1998; Liotta & Ranalli, 1999), indicating the difficulty of constraining this parameter. For the purpose of our simplified model, we have assumed a value of 10⁻¹⁴ s⁻¹, located in the middle of the range used for calculations by Liotta & Ranalli (1999).

Table 2. Rheological parameters related to the lithological scheme of Figure7.

Lithology	Rheological analog	Ac [MPa^−n s^−1]	n	E [J mol^−1]	Reference
Argillites	Wet quartzite	3.2 x 10^−4	2.3	1.54 x 10^5	Ranalli (1995)
Evaporites	Anhydrite	31.6	2	1.52 x 10^5	Muller & Briegel (1978)
Quartzites	Wet quartzite	3.2 x 10^−4	2.3	1.54 x 10^5	Ranalli (1995)
Phyllites	Biotite Schist	1.3 x 10^−67	31	0.98 x 10^5	Shea & Kronenberg (1992)
Gneiss	Wet granite	2.0 x 10^−4	1.9	1.37 x 10^5	Ranalli (1995)

Figure 7. Schematic lithological section of the Larderello geothermal area adopted for the rheological model: 1) Argillite; 2) Evaporite; 3) Quartzite; 4) Phyllite; 5) Gneiss.

5. Results of the thermal model

Size, depth and age of emplacement of the heat source, responsible for the regional anomaly of the surface heat flow (150–200 mWm⁻²) in the Larderello area, have been calibrated using a trial and error approach. Theoretical estimates have been compared with values of the surface heat flow observed along the NW-SE profile reported in Figure 4. After varying shape, size and depth of emplacement of the magmatic body, the preferred fitting was obtained with the following geometry of the heat source: a magmatic body having a flattened hexagonal shape in cross section; a thickness ranging from 8 km (at the centre) down to 4 km (at the border); a horizontal axis length of 14 km; a depth of the roof varying from a minimum of 3 km (at the centre) down to 5 km (at the border) (Figure 5). The regional heat flow anomaly is reproduced, satisfactorily, when 300 (± 100) kyr are elapsed since the magma emplacement (Figure 8). Instead, the peak value of the heat flow estimated by the model at the center of the regional anomaly (c. 300 mWm⁻²) turns out to be decidedly lower than the measured values, exceeding 500 mWm⁻² (see Figure 4) and locally reaching the peak of c. 1000 mWm⁻² (e.g., Bellani et al., 2004). This confirms that observed peak values above the pluton center reflect the dominant contribution of the convective heat transport (Norton & Knight, 1977). For this reason, our conductive model aims to reproduce values of the regional heat flow anomaly lower than 200 mWm⁻², related to zones close to the margins of the magmatic body.

Once constrained shape, size and depth of the pluton, the thermal evolution has been simulated from the surface down to 15 km, with several snapshots (Figure 9) during a magma cooling history of 1.4 Myr. Melt-present conditions are assumed for the zones with temperatures exceeding 650°C outlined by the purple contour line (Figure 9). In addition, values of the surface heat flow are provided at top of the contoured thermal sections.

Figure 9. Temperature (°C) cross-sections at various time steps after the magma emplacement. The white dashed line represents the boundary of the magmatic body. The purple line indicates the 650°C isotherm delimitating the melt-present zone. Model values of the surface heat flow (mWm⁻²) are also provided on the top of each cross-section (red arrows).

During the early thermal history (t < 200 kyr), the shallower part of the pluton is affected by a faster magma cooling with respect to zones close to the lateral margins. This reflects on the increase of the surface heat flow, quicker in the point located on the vertical axis of the pluton (max after 200 kyr) than in points located in correspondence of the pluton margins (max after 300 – 400 kyr). Closed contour lines are present up to 800 kyr. Then a slow depression and flattening of the isotherms reflects a relaxation towards a condition of thermal equilibrium. During the cooling history, the melt-present zone within the plutonic body progressively reduces and disappears at c. 1.2 Myr (Figure 9). In detail, this can be observed in the diagram of Figure 10a. Meanwhile, the upper front of the melt-present zone migrates from a minimum depth of 3 km to a depth of c. 11 km (Figure 10b). Interestingly, contoured thermal sections show that the zone below the pluton bottom may be affected by a significant heating, able to trigger partial melting in quartz-feldspars-rich rocks in water-saturated conditions. As previously stated, the snapshot at t = 300 kyr represents the best fit to the observed regional surface heat flow anomaly and thus to the current state of the Larderello area. The model predicts that the pluton roof has a temperature of 450 – 500°C and that the melt-present zone, reduced to 65% in volume, is located at a minimum depth of c. 5 km.

Figure 10. Melt reduction in the plutonic body as resulting in the thermal model. (a) Reduction and disappearance of the melt-present zone (vol %) during the cooling history of the magmatic body. (b) Migration of the upper front of the melt-present zone, from 3 km (depth of magma emplacement) to c. 11 km (at t = 1.2 Myr after magma emplacement).

We computed the evolution of the geotherm in the continental crust, along a vertical profile passing through the top vertex of the magmatic body (Figure 11a). It can be observed the thermal perturbation due to magma emplacement and the subsequent thermal relaxation, during the cooling history.

Figure 11. Diagrams illustrating the evolution of the transient geotherm and the variation of the surface heat flow reconstructed by the thermal model. (a) Evolution of the transient geotherm within the continental crust along a vertical profile passing through the top vertex of the magmatic body. The different colour lines indicate a specific time during the evolution of the transient geotherm. (b) Temporal variation of the surface heat flow calculated at the centre of the model anomaly.

We computed, also, the temporal variation of the surface heat flow at the center of the model anomaly (Figure 11b). It is important to note that the peak value is achieved 100–200 kyr after magma emplacement. The current situation is presumed to occur at 300 kyr after magma emplacement and corresponds to a phase of decreasing surface heat flow.

6. Results of the rheological model

Brittle and ductile strength profiles, calculated with the rheological model, have been used to obtain rheological sections at different time steps. Results are shown in Figure 12. In this figure, the BD¹ transition is represented by the black curve joining the points where the difference between brittle and ductile strength is equal to 0. It can be observed a significant rheological evolution with time and an evident migration of the BD¹ transition in proximity of the pluton. At t ≤ 300 kyr, the BD¹ transition coincides with the roof of the pluton. During the magma cooling, the ductile strength increases and the BD¹ transition deepens: at 800 kyr, it is within the magmatic body and at 1.4 Myr it tends to an equilibrium at a depth of c. 8.5 km. In Figure 13, the crustal strength profiles, along a vertical line crossing the centre of the plutonic body (x = 0 km), are represented at the same time steps as the rheological cross sections of Figure 12. The BD¹ transition, corresponding to the intersection between brittle and ductile strength profiles, is indicated by the blue and red lines, respectively. The decrease of the strength in the brittle domain, after c. 1.8 km, is due to the increase of λ with the depth.

Figure 12. Rheological cross-sections at various time steps after the magma emplacement. The colour scale denotes Log₁₀[strength, MPa]. The black line represents the location of the BD¹ transition. The white dashed line outlines the boundary of the magmatic body. The letters B and D indicate the brittle and ductile domains, respectively.

Figure 13. 1D strength profiles, along a vertical line crossing the centre of the plutonic body (x = 0 km), at different time steps after the magma emplacement. The blue line represents the brittle strength calculated by equation (7) and the red line represents the ductile strength calculated by equation (11). The dashed parts are related to probable melt-present conditions (T > 650°C).

Comparing the BD¹ transition profile with the K-horizon profile (Figure 14), an acceptable correspondence can be observed when 300 kyr are elapsed from the magma emplacement. In particular, this is evident in proximity of the roof of the pluton and on its right side. Instead, by choosing lower λ values, ranging between 0.80 and 0.95, it is observed a better correspondence on the left side of the pluton. Another possible influence on the BD¹ transition is linked to the strain rate value. Varying the exponent of the strain rate in the order of ± 1, with respect to the adopted value (10⁻¹⁴ s⁻¹), the BD¹ transition’s depth is approximately unchanged in correspondence of the pluton roof and it changes laterally, of c. ± 1.5 km.

Also, comparing the thermal section at 300 kyr with the K-horizon and BD¹ transition profiles (Figure 15), it can be observed that their depths, near the pluton roof, correspond broadly with the 450–500°C isotherms. Instead, on the pluton sides, the BD¹ transition profile coincides with the 400°C isotherm. Finally, the K-horizon profile follows the 300–350°C isotherms on the left side and the 350–400°C isotherms on the right side of the magmatic body.

7. Discussion

The difficulty to understand the relationship between the magmatic bodies and geothermal systems has been underlined in several papers (Cathles et al., 1997; Erkan, Blackwell, & Leidig, 2005; Stimac et al., 2001). This can be particularly critical where evidence of volcanic activity is poor or absent, as in the Larderello area (Erkan et al., 2005; Mongelli et al., 1998). In our case, the magmatic heat source is not clearly defined: size, age and depth of emplacement of the magmatic body have not been constrained yet. The thermal model presented in this paper, despite the simplifications and the shortcomings to represent completely a geothermal system, may be useful to cover these gaps of knowledge. The choice of a purely conductive model can be allowed if the objective is the reproduction of the regional thermal anomaly due to the emplacement of a magmatic body in the upper crust. This approach has been adopted by Dalrymple et al. (1999), for the case of the Geysers geothermal field.

Another simplification in the model is the assumption of the instantaneous magma emplacement, in contrast with the growing evidence that emplacement of medium to large magmatic bodies occurs incrementally (e.g. Annen, 2011). However, for the case of the Larderello area, geometric and chronological data on the plutonic body are too scanty to define realistically size and frequency of the magma increments. For this reason, an instantaneous emplacement represents a forced choice, that still can provide significant results for end-member conditions. The thermal evolution for a single intrusion depends, on size and shape of the pluton, magma and crust initial temperatures, depth of magma emplacement and thermal conductivity of the host rocks (Annen, 2011; Stimac et al., 2001). In our case, varying the size and depth of emplacement, all other factors being equal, the most important key parameter influencing the surface heat flow is the depth of emplacement of the magma body. In fact, the peak value of the regional thermal anomaly cannot be reproduced if the apex of the pluton roof is located at a depth exceeding 3 km. Stimac et al. (2001), reached a similar conclusion in the study of the Clear Lake magmatic hydrothermal system (California, USA), finding a depth range of the pluton of 3–4 km, in order to reproduce the observed surface heat flow anomaly. Interestingly, the estimated depth of emplacement of the magmatic body, from 3 (top vertex) to 5 km (lateral vertices) is in agreement both with geological (Bertini et al., 2006; Dini et al., 2005), geophysical (Fiordelisi et al., 1995; Gianelli et al., 1997) and petrological (Gianelli, 1998) data available in literature. Regarding the size, the estimated horizontal axis length of the pluton (14 km) is in agreement both with the width proposed by Mongelli et al. (1998), on the basis of the magnetotelluric data (Fiordelisi et al., 1995), and compatible with the dimension of the pluton deduced from the geological cross-sections published by Bertini et al. (2006). The estimated age of emplacement of the magmatic body (300 kyr) is in agreement with the model estimates of Mongelli et al. (1998) and Gola et al. (2017) and turns out to be comparable with the age estimated for the volcanic activity of the neighbouring Mt. Amiata (Dini et al., 2005). However, it is decidedly lower than geochronologic dates (> 1 Myr) available for some intrusive rocks cored in exploratory wells in the Larderello geothermal field (see Table 1). The discrepancy between geochronological dates and numerical model results was found by several authors (e.g., Cathles et al., 1997), for different geothermal areas in the world. It represents a critical point for a complete comprehension of the relationship between geothermal systems and heat sources. The problem cannot be solved by considering the effect of the convective heat transfer. In fact, in this case, a faster cooling of the magma body would be obtained and thus an even younger age of emplacement would be estimated. The same effect would have the implementation of a 3D model, implying a faster heat loss of the magmatic body than a 2D model. A chance to explain the age mismatch is provided by models based on the incremental emplacement of the magmatic body (Annen, 2011). By this way, it would be possible to reproduce a more persistent thermal anomaly. However, to avoid a drastic reduction of the simulated surface heat flow, the interval between single magmatic pulses should be contained below a critical value (e.g., 100 kyr; see later). For the same reason, an upward accretion of the magmatic body is preferred with a respect to a downward accretion. On the basis of these considerations, an incremental growth of the plutonic body was reproduced by a numerical model (see Langone, Caggianelli, Festa, & Prosser, 2014 for details). The adopted approach consisted in finding the maximum time interval between single magma pulses still compatible with the observed values of the regional surface heat flow anomaly. Indeed, if the time interval between single magma pulses is much greater than the solidification time of a single pulse, i.e. in case of slow pluton growth (Annen, 2011), the local thermal perturbation generates a negligible effect on the surface heat flow. In the absence of chronological and geometric data, the number of the magma increments, has been constrained by the need to preserve the melt-present conditions during the pluton growth and by the final size of the modelled magmatic body. Then, considering a pluton growth occurring in ten magma increments, it was obtained a maximum interval between single magma pulses of c. 100 kyr, that implies a pluton growth occurring in 1 Myr. By this way, the observed value of the surface heat flow, in correspondence of the pluton margins (average value of 170 mWm⁻²), can be still reproduced (Figure 16). By increasing the time interval to 200 and 300 kyr and a pluton accretion completed in 2 Myr and 3 Myr, respectively, the peak of the surface heat flow progressively decreases being well below the observed reference value (see Figure 16). Therefore, by this approach the maximum age estimated for the start of the emplacement of the Larderello granite is of 1 ± 0.3 Myr, lower than most of the available geochronological dates. The only comparable value of 1.23 ± 0.13 Ma (see Table 1) has been obtained by Villa, Ruggieri, & Puxeddu (2001). On the basis of our results, the other magmatic rocks that provided ages older than 2 Myr cannot be considered responsible for the current thermal anomaly of the Larderello area. It is worth noting that, even considering an incremental growth, the observed heat flow is reproduced only in the final stages of pluton accretion. Then the results of both models (single and incremental emplacement) suggest recent last injections of magma, in the Larderello area. This result has important implications for understanding the geothermal field, because it demonstrates that although there have been various phases between 3.8 and 0.3 Ma in the area, we can consider that there is a relative continuity in this magmatic activity.

Figure 16. Variation of the surface heat flow with time, in correspondence of the pluton margin, according to the incremental growth model of the magmatic body, in case of a pluton construction in 1, 2 and 3 Myr. The horizontal dashed line indicates the corresponding value of the observed surface heat flow. Magma increments have volume increasing from c. 70 km³ up to a maximum of 250 km³ and then decreasing back to 70 km^3.

If the pluton growth was accomplished in an interval of c. 1 Myr, an estimate of the rate of pluton construction can be made. Firstly, a volume of c. 2000 km³ can be calculated for the Larderello plutonic body, by considering the 2D section of the magmatic body assumed in the model and its probable extent in NE-SW direction (see Figure 4). It transpires that the pluton construction rate was of c. 2. x 10⁻³ km³yr⁻¹. This value is compatible with data on duration and pluton construction rates compiled by De Saint Blanquat et al. (2011), ranging between 10⁻¹ and 10⁻⁴ km³yr⁻¹. Nonetheless, further studies, dedicated to model more complex mechanisms of emplacement and growth of the Larderello granite, are needed for a full comprehension of the connection between the regional thermal anomaly and the magmatic heat source.

Regarding the debated interpretation of the K-horizon, despite the simplified assumptions involved in the estimation of the 2D strength profiles, the results show an acceptable correspondence between the BD¹ transition and the K-horizon profile when 300 kyr are elapsed from the magma emplacement. This is compatible with the hypothesis that the K-horizon represents the top of the BD¹ transition (e.g., Liotta & Ranalli, 1999). Comparing the thermal cross-section at 300 kyr with the K-horizon profile it results that it occurs at a temperature of c. 400 ± 100°C, in accordance with boreholes data showing that the seismic reflector follows approximately the 450 ± 50°C isotherm (e.g., Liotta & Ranalli, 1999). Moreover, from the evolution of 2D strength profiles of the upper crust, it emerges that the BD¹ transition is time-dependent and that its pronounced culmination progressively fades out during magma cooling. Throughout this migration there is a time interval (up to 300 kyr) when it coincides with the pluton roof, conforming to the interpretation of the K-horizon provided by Bertini et al. (2006) and more recently by Saccorotti et al. (2014), on the basis of local earthquake tomography.

Therefore, from our results it emerges that the interpretations proposed of the K-horizon, i.e. (i) the top of the BD¹ transition (e.g., Liotta & Ranalli, 1999) vs. (ii) the roof of the pluton (e.g. Bertini et al., 2006), are not necessarily in disagreement.

The presence of extensional shear zones in correspondence with the regional BD¹ transition may have played an important role for the rising of granitic magma in the crust. It is worth noting that the link between extensional regime and syn- to late-tectonic plutonism has been suggested, also, for the emplacement of the Porto Azzurro pluton in the near Elba Island (Italy) (Smith, Holdsworth, & Collettini, 2011), which is considered as an analogue of the current Larderello geothermal field (e.g., Dini, Mazzarini, Musumeci, & Rocchi, 2008). The ascent of further magma above the regional BD¹ transition, may be favoured by nontectonic processes, accommodating the space required for the pluton emplacement, like roof lifting (Acocella & Rossetti, 2002), which we retain could explain the upper shape of the modelized magmatic body. This scenario may be compatible with our results, indicating that the roof of the pluton coincides with the current position of the local BD¹ transition. Nevertheless, specific studies are necessary for understanding the relationship between deformation and pluton emplacement in geothermal areas.

8. Concluding remarks

The proposed thermo-rheological model allowed us to constrain size, age and depth of emplacement of the magmatic heat source responsible for the regional heat flow anomaly of the Larderello area. In addition, rheological calculations provided information useful for the interpretation of the K-horizon, a seismic reflector of controversial interpretation.

The main results of the present study are summarized below:

• The regional surface heat flow anomaly can be ascribed to the emplacement of an approximately biconvex magmatic body having a maximum thickness ranging from 8 km (at the centre) down to 4 km (at the border), a horizontal axis length of 14 km and a depth of the roof varying from a minimum of 3 km (at the centre) down to 5 km (at the border).

• The anomaly is reproduced satisfactorily when 300 kyr are elapsed from the magma emplacement. At this time, the pluton roof has a temperature of 450–500°C, the melt- present zone reduced to 65% of the original volume and the upper magmatic front retreated to a minimum depth of c. 5 km.

• The estimated age of 300 kyr is decidedly lower than the current geochronologic dates (> 1 Myr) of the Larderello intrusive rocks (Dini et al., 2005). Also considering an incremental emplacement of the magmatic body the age mismatch cannot be explained. A maximum of 1 ± 0.3 Ma is required for the start of the pluton growth. On this basis, it is sustained that magmatic rocks below the Larderello area, older than 1.3 Ma cannot be responsible for the current regional heat flow anomaly.

• The pluton crystallization is completed in 1.2 Myr when the surface heat flow anomaly becomes negligible.

• The X-Z mechanical strength profiles allowed to reconstruct the drifting of the BD¹ transition during the cooling of the magmatic body.

• The K-horizon profile is compatible with both the upper contact of the plutonic body and the BD¹ transition, when 300 kyr are elapsed from the magma emplacement.

Finally, the performed study, although applied to a specific case, may provide interesting results for all the geothermal areas related to magma emplacement in the upper continental crust.

Acknowledgments

We are grateful and in debt with Giorgio Ranalli for showing the way to build up a thermo-rheological model. Discussions with Domenico Liotta were useful to set up the rheological model and to refine the interpretation of our results. Comments and suggestions by the editor Erdin Bozkurt and an anonymous reviewer helped to improve significantly the former version of the manuscript.

Disclosure statement

No potential conflict of interest was reported by the authors.