Controls of Radiogenic Heat and Moho Geometry on the Thermal Setting of the Marche Region (Central Italy): An Analytical 3D Geothermal Model

Abstract

Using published cross-sections and a series of geological constraints, a 3D geological model of an important area of the Adriatic sector of peninsular Italy—i.e., the Marche region—was developed. Then, an analytical procedure, taking into account the heat rising from the mantle and the radiogenic heat produced by the crust, was applied on the pre-built structural model, in order to obtain the 3D geothermal setting of the entire region. The results highlighted the key role played by the Moho geometry, particularly as a step of ~10 km occurs between the Adriatic Moho of the subducting plate to the west and the new Tyrrhenian Moho characterizing the back-arc area to the west. The comparison between our results and available borehole data suggests a good fit between the applied analytical methodology and published datasets. A visible anomaly is located at a specific site (i.e., the coastal town of Senigallia), where it may be envisaged that fluid circulation produced a local surface heat flow increase; this makes the Senigallia area a promising feature for the possible exploitation of geothermal systems.

1. Introduction

The high geothermal energy potential is well known in central Italy, especially in the Tuscany region, where geothermal exploration started during the 1990s in Larderello and Monte Amiata [1,2,3,4,5,6]. At present, closer attention has been given to the geothermal energy potential of the neighboring Marche region, involving the investigation of both fluid circulation [7] and surface heat flux [8,9,10]. Recent studies focused on the central Apennines regional setting, which is characterized by an increasing surface heat flow moving from the mountain chain to the Adriatic off-shore to the east, besides that occurring toward the Tuscany region to the west [8,9,10].

When assessing a potential site for geothermal exploitation, a comprehensive picture of crustal setting and temperature profile is fundamental for any site-specific, appropriate field development. Today, indispensable geothermal suitability assessments are often carried out using a series of invasive inspections involving legal permission and high cost (e.g., well drilling) [11].

In this study, we produced the first 3D geological model of the crust in the Marche region based on Moho depth provided by Grad et al. [12], the seismic profile CROP-03 interpreted by Barchi et al. [13] and Santini et al. [9], a series of published balanced geological cross-sections [14,15,16], the geological map by Conti et al. [17] and a 10 m cell size digital elevation model (DEM) [18]. Subsequently, taking into account the heat rising from the mantle and the radiogenic heat produced by radioactive elements in the crust, an analytical methodology was used in order to produce a crustal 3D thermal model of the entire sector based on the pre-built 3D geological model. Finally, we compared our results with those obtained by previous studies to improve the description of thermal implications. This work provides a general framework of temperature distributions in the studied area, and emphasizes the key role played by the Moho geometry in the thermal structure of the central Apennines.

2. Geological Background

The area of study (Figure 1) is located in the Apennine fold and thrust belt. This area represents the most recent expression of the geodynamic process that produced the Mediterranean basins after the collision between Europe and Africa, which started during the late Cretaceous period and was followed by Oligocene and Neogene periods’ slab-retreat events [15,17,19,20,21,22].

The tectonic evolution of the entire area is strongly debated due to the complexity of the geodynamic processes that generated the mountain chain (e.g., [9,21,23,24]). Several authors have proposed models based on stratigraphic and structural analysis of transects oriented along the flow lines of relative motion [25,26], whereas others have followed an approach based on the laws of plate kinematics [22,27]; however, all of these studies found a Moho step in correspondence with the Apennines chain, caused by a discontinuity between the subducting Adriatic plate and the Tyrrhenian back-arc domain to the west. The Moho step (~10 km) was also observed by recent geophysical studies of Italian and European Moho geometry [12,21,28,29,30,31,32,33].

In this study area, the deepest units are the Paleozoic crystalline basement and the overlaying Permo-Triassic continental siliciclastic strata (i.e., the Verrucano Group). Neither of them crop out in the area, but their presence was confirmed by seismic studies and deep wells [28,34]. The overlying sedimentary cover is several thousand meters thick. It is composed of Triassic evaporites (Burano Anhydrites Fm.) followed by the ca. 1000 m thick, pre-orogenic and dominantly carbonate, Triassic–Eocene portion of the Umbria–Marche succession (including the Calcare Massiccio Fm. at the base, up to the Scaglia Rossa Fm. at the top). This is overlain by a 2000–2500 m thick succession of Late Eocene to Plio-Pleistocene, dominantly terrigenous marine deposits (including the Scaglia Cinerea Fm. up to the Sabbie Gialle Fm.), followed by late Quaternary continental deposits (Figure 1) [35,36,37,38].

Overlying the Triassic evaporites, the Late Triassic–Lower Jurassic platform carbonates of the Calcare Massiccio Fm. reach 800 m of thickness. During the Jurassic extensional phase, the carbonate platform was dissected by normal faults, producing structural depressions and structural heights [39]. On this uneven bathymetry, pelagic sediments of variable thickness were deposited forming ‘complete’, ‘condensed’ and ‘composite’ successions [40,41]. Basin deposition of limestone and marl continued from Late Jurassic to Middle Miocene times, progressively increasing the contribution of argillaceous material in younger formations. The closure of the Mesozoic Tethys Ocean and the collision of European (Corsica–Sardinia block) and African (Adria block) continental margins, resulted in collisional orogenesis since the Oligocene [42]. The resulting fold-and-thrust belt consists of ENE-verging folds and WSW-dipping thrusts [14,15,16,43,44,45,46,47,48,49]. Coeval deposition of siliciclastic material occurred in the evolving foreland basin system (including, in the Marche region, the Miocene–Arenacea, Camerino and Laga basins [50], as well as the Plio-Pleistocene Argille Azzurre basin (stratigraphically overlying Messinian evaporites). Late Quaternary continental deposits unconformably overly the marine succession.

The stress field acting in the axial zone of the central Apennines since the Middle Pleistocene produced normal fault systems responsible for moderate to high seismicity [51,52] in the inner zone of the Marche region [16,53,54,55,56,57,58,59,60,61,62,63,64,65,66], whereas active horizontal compression characterizes the coastal sector and offshore areas [67,68].

3. Materials and Methods

In order to produce a 3D geothermal model of the Marche region, we used an analytical procedure based on a 3D geological model and calculated geotherms on a series of pseudo-wells located within the 3D geological model. The results for each pseudo-well were then interpolated to produce a 3D model of the entire thermal structure of the region.

3.1. 3D Geological Model

Over the years, many authors produced geological cross-sections perpendicular to the strike of the mountain chain based on available seismic profiles and geological field observations [7,13,14,15,16,21,49,69]. Furthermore, seismic tomographies were performed to investigate the thickness of the sedimentary cover and the Moho depth [12,21,29,31,32,33].

We produced a three-layer 3D geological model of the study area based on the Moho geometry provided by Grad et al. [12], the inner zone of the seismic profile CROP-03 [13] interpreted by Mazzoli et al. [14] (Figure 2), a series of published balanced geological cross-sections [14,15,16] (Figure 3), the geological map by Conti et al. [17] (Figure 1), and a 10 m cell size digital elevation model (DEM) [18].

The three-layer 3D geological model (Figure 4) was built using Blender [70], a free and open-source 3D computer graphics software used to show complex 3D geometries in many fields of study [71,72], including geology and geophysics [73,74,75].

3.2. Geothermal Model: Constraints and Assumptions

The computational procedure took into account the heat rising from the mantle and the radiogenic heat due to radioactive elements located in the crust. The latter is divided, in this case, into three layers (Figure 4). To perform the geothermal model, we considered one hundred and sixty-five pseudo-wells located in the nodes of a ca. 10 × 10 km sampling grid traced on the study area, in addition to other six pseudo-wells placed on CROP-03 (Figure 5).

Each pseudo-well included data on the heat source, Moho depth, and thickness of each layer considering:

1

Altitude: we took into account the topography using a 10 m cell size digital elevation model (DEM) [18];

2

Variable thickness of the first layer $\left(h_{CS}\right)$ represented the terrigenous marine succession (Late Eocene to Plio-Pleistocene) and late Quaternary continental deposits. It was based on a series of published balanced geological cross-sections [14,15,16] (Figure 3), and the geological map by Conti et al. [17] (Figure 1). We considered a constant heat production rate of the Siliciclastic succession $H_{CS}=1.05 \mathsf{\mu }\mathrm{W}{\mathrm{m}}^{-3}$ (with a range of $0.9–1.3 \mathsf{\mu }\mathrm{W}{\mathrm{m}}^{-3}$) [10] and thermal conductivity $k_{CS}=2.1 \mathrm{W}{\mathrm{m}}^{-1} °{\mathrm{C}}^{-1}$ (with a range of $2.0–2.2 \mathrm{W}{\mathrm{m}}^{-1} °{\mathrm{C}}^{-1}$) [8];

3

The variable thickness of the second layer $\left(h_{CC}\right)$, representing the Umbria–Marche calcareous–marly succession from Triassic evaporites (Burano Anhydrites Formation) to the Scaglia Rossa Fm., was based on a series of published balanced geological cross-sections [14,15,16] (Figure 3), and the geological map by Conti et al. [17] (Figure 1). We considered the constant heat production rate of the Umbria–Marche calcareous–marly succession, $H_{CC}=0.45 \mathsf{\mu }\mathrm{W}{\mathrm{m}}^{-3}$ (with a range of $0.3–0.6 \mathsf{\mu }\mathrm{W}{\mathrm{m}}^{-3}$) [10], and thermal conductivity, $k_{CC}=2.4 \mathrm{W}{\mathrm{m}}^{-1} °{\mathrm{C}}^{-1}$ (with a range of $2.3–2.5 \mathrm{W}{\mathrm{m}}^{-1} °{\mathrm{C}}^{-1}$) [8];

4

Variable thickness of the third layer $\left(h_B\right)$ represented the basement. We considered a thickness of the basement from the Moho (based on Grad et al. [12]) to the top basement constrained by the CROP-03 seismic profile [13] interpreted by Mazzoli et al. [14] (Figure 2) and a series of published balanced geological cross-sections [14,15,16] (Figure 3). We considered a basement thermal conductivity of $k_B=2.7 \mathrm{W}{\mathrm{m}}^{-1} °{\mathrm{C}}^{-1}$ (with a range of $2.6-2.8 \mathrm{W}{\mathrm{m}}^{-1}°{\mathrm{C}}^{-1}$) [8] and a heat production rate of the basement of $H_B=3.15 \mathsf{\mu }\mathrm{W}{\mathrm{m}}^{-3}\left(\mathrm{with}\mathrm{a}\mathrm{range}\mathrm{of}3.0-3.3 \mathsf{\mu }\mathrm{W}{\mathrm{m}}^{-3}\right)$ [76,77], exponentially decreasing with depth. The radiogenic sources’ intensity, for any thickness $h_B$, diminished downward with a logarithmic decrement $D$, which had the dimension of depth and a characteristic value of 8 km in the region [76,78];

5

Variable heat flux at the Moho $\left(Q_m\right)$ from $20 \mathrm{m}\mathrm{W}{\mathrm{m}}^{-2}$ to the east and $40 \mathrm{m}\mathrm{W}{\mathrm{m}}^{-2}$ to the west [21,79,80].

3.3. Analytical Procedure

The entire thermal structure is affected by two sources of heat: (i) the heat flowing upward from the Moho, coming from the mantle $\left(Q_m\right)$, and (ii) the heat production rates of the cover ($H_{CS}$ and $H_{CC}$) and the basement $\left(H_B\right)$, due to the presence of radioactive elements in the crust. To obtain the surface heat flux, we used the analytical procedure by Dragoni et al. [76] and further applied by Santini et al. [9,78] and Basilici et al. [74,81], modified specifically for the area of study. In the computation of the surface heat flow, in addition to contributions of $Q_m$, we added radiogenic heat, obtaining the following equation used to calculate the surface heat flow $\left(Q_s\right)$:

$$Q_S=Q_m+H_{CS}{h}_{CS}+H_{CC}{h}_{CC}+H_{B}D\left(1-e^{-\frac{\left(h-h_C\right)}{D}}\right)$$

(1)

where $h$ and $h_C$ are the thicknesses of the crustal structure and of the whole cover, respectively, with $h_C=h_{CS}+h_{CC}$.

To compute the temperature trend with the depth, we adopted the analytical procedure summarized in

Table 1. Equations adopted for the analytical procedure.

	0 ≤ z ≤ hCS	hCS < z ≤ hC	hC < z ≤ h
$T\left(z\right)$	$T_s+T_{CSH}\left(z\right)$	$T_s+T_{CSH}\left(h_{CS}\right)+T_{CCH}\left(z\right)$	$T_s+T_{CSH}\left(h_{CS}\right)+T_{CCH}\left(h_C\right)+T_{BH}\left(z\right)$
$T_{CSH}\left(z\right)$	$\frac{Q_S}{k_{CS}}z-\frac{H_{CS}}{2k_{CS}}{z}^2$
$T_{CCH}\left(z\right)$		$\frac{Q_S-H_{CS}{h}_{CS}}{k_{CC}}\left(z-h_{CS}\right)-\frac{H_{CC}}{2k_{CC}}{\left(z-h_{CS}\right)}^2$
$T_{BH}\left(z\right)$			$\frac{Q_m }{k_B}\left(z-h_C\right)+\frac{H_B{D}^2}{k_B}\left(1-e^{-\frac{\left(z-h_C\right)}{D}}\right)$

, while the parameters considered are listed in

Table 2. Analytical procedure parameters.

$h_{CC}$	calcareous cover thickness
$h_C$	whole cover thickness
$h_B$	basement thickness
$h$	whole crustal thickness
$D$	depth scale
$H_{CS}$	siliciclastic cover radioactivity produced heat
$H_{CC}$	calcareous cover radioactivity produced heat
$H_B$	basement radioactivity produced heat
$k_{CS}$	siliciclastic cover thermal conductivity
$k_{CC}$	calcareous cover thermal conductivity
$k_B$	basement thermal conductivity
$Q_s$	surface heat flow density
$Q_m$	mantle heat flow density
$T_s$	surface temperature
$T_{CSH}$	siliciclastic cover heat temperature
$T_{CCH}$	calcareous cover heat temperature
$T_{BH}$	basement heat temperature
$h_{CS}$	siliciclastic cover thickness

(with symbols and related definitions).

The temperature in the first layer $\left(h_{CS}\right)$ is due to the heat flow at the top of calcareous cover and to the contribution of the radiogenic source of this layer, adding to these, the surface temperature ($T_S$)

$$T\left(z\right)=T_S+\frac{Q_S}{k_{CS}}z-\frac{H_{CS}}{2k_{CS}}{z}^2,0\le z\le h_{CS}$$

(2)

for simplicity, defining the depth dependent terms as $T_{CSH}\left(z\right)$.

In the second layer $\left(h_{CC}\right)$, the geotherm initial value is the temperature at the bottom of the first layer $\left(T_S+T_{CSH}\left(h_{CS}\right)\right)$, to which is added the temperature variation due to the radiogenic source present in the calcareous cover and to the heat flow at its bottom

$$T\left(z\right)=T_S+T_{CSH}\left(h_{CS}\right)+\frac{Q_S-H_{CS}{h}_{CS}}{k_{CC}}\left(z-h_{CS}\right)-\frac{H_{CC}}{2k_{CC}}{\left(z-h_{CS}\right)}^2,h_{CS}\le z\le h_C$$

(3)

defining the depth dependent terms as $T_{CCH}\left(z\right)$.

Finally, the temperature in the basement depends on the mantle heat flow and on the radiogenic source present in it, adding the temperature calculated at the top of the layer ($T_s+T_{CSH}\left(h_{CS}\right)+T_{CCH}\left(h_C\right))$

$$T\left(z\right)=T_S+T_{CSH}\left(h_{CS}\right)+T_{CCH}\left(h_C\right)+\frac{Q_m }{k_B}\left(z-h_C\right)+\frac{H_B{D}^2}{k_B}\left(1- e^{-\frac{\left(z-h_C\right)}{D}}\right),h_C\le z\le h$$

(4)

4. Results

Using the analytical procedure described in Section 3.3, we calculated surface heat flow, geotherms, and isotherms for each pseudo-well (Figure 5) located in the 3D geological model (Figure 4). The output model of the surface heat flow computation is shown in Figure 6. It was obtained by interpolating the $Q_S$ value obtained from each pseudo-well with a fourth-order polynomial equation.

The variable thickness of the three layers of the 3D geological model produced a variable trend in the surface heat flow, which appears to be strongly influenced by the Moho geometry. Surface heat flow varies from $46–50 \mathrm{m}\mathrm{W}{\mathrm{m}}^{-2}$ in the central zone, to $60–67 \mathrm{m}\mathrm{W}{\mathrm{m}}^{-2}$ in the inner zone (Arezzo–Perugia area). This variation appears to be mainly associated with the Moho step, which rises the top of the mantle by ~10 km in the SW sector (Figure 4a). In addition, the Rimini–Cesena area and the zone located between Ancona and Macerata, show an increase in surface heat flow of ~5 $\mathrm{m}\mathrm{W}{\mathrm{m}}^{-2}$, probably due to the thickening of the Siliciclastic succession (Figure 4b,c). Furthermore, considering the range of values of the radiogenic heat sources for the three layers, we estimated the changes in the computed surface heat flow. In particular, we obtained a maximum increase of about 6% if the upper limits of the range are considered, whereas the decrease in surface heat flow did not exceed 5%, taking into account the lower limits of the range.

To study the trend in the temperature and depth along the CROP-03 line, we calculated the geotherms for the six pseudo-wells located along the profile, and then interpolated them using second-order polynomial equations representing the 50°, 100°, 150° and 200° isotherms. Figure 7 and Figure 8 show the geotherms calculated for the pseudo-wells CR1 to CR6 (Figure 5) and the relative isotherms along the CROP-03 (Figure 2), respectively. In particular, pseudo-well CR2 showed an overall higher temperature profile with respect to CR1 (Figure 7). This was due to the higher altitude (of 1000 m a.s.l.) of CR1 with respect to the other pseudo-wells (located 0–400 m a.s.l.). Nevertheless, the general trend in the calculated geotherms and interpolated isotherms shows an increase in the temperature gradient going from the coastal area to the inner zone.

Furthermore, in order to analyze the trend in the isotherms for the whole study area, we calculated the geotherm for each pseudo-well, and then interpolating by second-order polynomial equations, we obtained the contour maps representing the 50 °C, 100 °C, 150 °C and 200 °C isotherms (Figure 9). The contour map for the isotherm at 50 °C shows a general, northward deepening trend characterized by an area of shallow depth (<1.1 km), south of Camerino. On the other hand, the isotherm 100 °C map shows a dominantly north-eastward deepening trend, with its shallowest depths in the westernmost area (close to Arezzo). The contour maps of 150 °C and 200 °C isotherms show similar trends characterized by greater depths in a wider area—defined as a ‘thermal depression’—extending from Pesaro to Camerino and just west of Macerata, whereas the shallowest depths occur in the westernmost area (close to Arezzo and Perugia).

5. Discussion

The thermal structure of the area is characterized by a trend that highlights a general increase in temperature from the coastal area to the inner zone. This is visible in the model of the surface heat flow (Figure 6), along the CROP-03 profile (Figure 7 and Figure 8) and in the temperature contour map (Figure 9). The analytical procedure indicates that the geometry of the Moho produces a high impact on the entire thermal setting of the Marche region, showing high heat flow values where the Moho is shallower (Arezzo–Perugia–Foligno area) and a sharp decrease in heat flow values, corresponding to the ~10 km Moho step. The presence of radioactive elements in the crust, instead, produces a homogeneous increase in temperature in the entire area, especially where the siliciclastic layer increases in thickness. However, the difference in the heat production rates between the calcareous layer $(H_{CC}=0.45 \mathsf{\mu }\mathrm{W}{\mathrm{m}}^{-3}$) and the siliciclastic layer $\left(H_{CS}=1.05 \mathrm{m}\mathrm{W}{\mathrm{m}}^{-3}\right)$ is substantial. A difference of $0.60 \mathrm{m}\mathrm{W}{\mathrm{m}}^{-3}$ in the heat production rate is also responsible for the ‘thermal depression’ located between Pesaro, Camerino and Macerata, where the siliciclastic cover is almost absent. On the other hand, the rise of the 50 °C isotherm south of Camerino could be linked with the circulation and ascent towards the surface of fluids in correspondence with the most important calcareous reliefs.

To obtain a better understanding of the thermal structure of the study area, our results were compared with previous studies [8,10,82]. Pauselli et al. [10] produced a surface heat flow map of central Italy constrained by available thermal and petrophysical logs from wells drilled for geothermal and hydrocarbon exploration purposes. In Figure 10, we compared our surface heat flow map with the results of Pauselli et al. [10].

The general trend in the surface heat flow calculated by the improved analytical procedure is comparable with that obtained by Pauselli et al. [10]. The Arezzo–Perugia area represents the highest value of $Q_S$ in the entire region (~$70 \mathrm{m}\mathrm{W}{\mathrm{m}}^{-2})$ in both our study (Figure 6) and in Pauselli et al. [10] (Figure 10). Slightly more elevated values in the map created by Pauselli et al. [10], with respect to our result, suggest a possible source of heat that is in addition to the radiogenic heat produced by the crust and the heat flowing up from the mantle $\left(Q_m\right)$, which are included in the analytical procedure. In the Senigallia area, a $Q_S$discrepancy value is clear. This means that in Senigallia, the surface heat flow measured from the wells (Figure 10) is due to another heat source. Chicco et al. [7] highlighted the presence of deep, high-angle faults that can certainly influence fluid circulation by means of water ascent from depth along fault damage zones. This fact makes the Senigallia area a promising feature for possible geothermal systems. Other $Q_S$ discrepancies occur in the Rimini area, although they are less prominent. The rest of the map does not show $Q_S$ discrepancies, thus confirming that most of the measured surface heat flow is due to the radiogenic heat and the heat flowing up from the mantle, which, in turn, depends on Moho geometry. This fact is visible in the Camerino–Pesaro area, where the absence of the siliciclastic cover results in a lower $Q_S$ value in both the analytical model (our results, Figure 6) and measured values (Pauselli et al. [10] result, Figure 10).

Della Vedova et al. [8] computed isotherms along the CROP-03 profile using boreholes within a 10 km distance from the section trace and lithology-dependent thermal conductivity values. The result of the work is compared with our results in Figure 11.

The general trend pointed out by the analytical procedure (blue isotherms in Figure 11) is consistent with the results by Della Vedova et al. [8]. The temperature values discrepancy is probably due to a different mathematical interpolation used. Obviously, the T value discrepancy increases with depth. To improve the description of our results, in Figure 12, we compared three well temperature profiles published by Pauselli et al. [10], Della Vedova et al. [8] and Verdoja et al. [81] with the closest pseudo-wells of our study.

The comparison between our pseudo-wells (orange lines in Figure 12) and published temperature profiles (green lines in Figure 12) shows a good fit between the analytical model and borehole measurements.

6. Conclusions

The results of this study provide the first 3D geological model of the crust in the Marche region, based on Moho depth by Grad et al. [12], the deep seismic reflection profile CROP-03 [13] interpreted by Mazzoli et al. [14], a series of published balanced geological cross-sections [14,15,16], the geological map by Conti et al. [17] and a 10 m cell size digital elevation model (DEM) [18]. An analytical methodology was elaborated to produce a 3D geothermal model of the crust taking into account the heat rising from the mantle and the radiogenic heat produced by radioactive elements. The model output showed an increase in temperature going from the coastal area to the inner zone. Comparing our outcomes with published results of previous works [8,10,81], we highlighted compatible results between our model and borehole measurements. The Moho geometry appeared to play a key role in the temperature model of the entire Marche region, producing a sharp change in temperature where a step of ~10 km occurs between the Adriatic foreland and Tyrrhenian back-arc domains. In the Senigallia and Rimini areas, a discrepancy of surface heat flow occurred between our $Q_S$ values and those provided by Pauselli et al. [10], probably due to circulating and rising fluids [7] not considered in the applied analytical procedure; this fact makes the Senigallia area a promising site to evaluate possible geothermal systems. In particular, besides the suitable condition for low-enthalpy plants for the heating/cooling of buildings, dedicated studies should be aimed at investigating the possibility of medium-enthalpy system plants.