Thermal Structure of the Northern Outer Albanides and Adjacent Adriatic Crustal Sector, and Implications for Geothermal Energy Systems

Using an analytical methodology taking into account heat flow density data, frictional heating, temperature variations due to the re-equilibrated conductive state after thrusting and geological constrains, we calculated surface heat flow, geotherms and isotherms along a balanced and restored regional geological cross-section. Our results highlight the impact of frictional heating produced by thrusts on the thermal structure of the study area, leading to a raising of the isotherms both in the inner Albanides to the E and in the Adriatic sector offshore. Minimum values of Qs in the surroundings of Tirana and the reconstructed 2D thermal structure suggest less favorable conditions for exploitation of geothermal energy, besides the direct use (Borehole Heat Exchanger-Geothermal Heat Pump systems). Nevertheless, the occurrence of the “Kruja geothermal zone”, partially overlapping this area and including hot spring manifestations, emphasize the structural control in driving hot fluids to the surface with respect to the regional thermal structure.


Introduction
Based on a new geological cross-section, an analytical methodology was implemented to produce a crustal thermal model that takes into account a series of geologically derived constraints and the temperature variation due to the re-equilibrated conductive state associated with thrusting, as well as heat flow density data and frictional heating. The resulting model obtained for the study area provides thermal constraints that are fundamental for the assessment and management of energy resources (hydrocarbons, geothermal energy) and also constitute a basis for tectonic and geodynamic modeling. In fact, knowing the thermal structure of the continental lithosphere is pivotal for the study of tectonic, geodynamic and surface processes, e.g., [1][2][3][4]. At present, a series of temperature profile measurements from selected boreholes are available for the entire country of Albania, resulting in a regional setting characterized by an increasing geothermal gradient from the west (Peri-Adriatic Depression) to the east (Internal Albanides), although its detailed pattern in the inner Albanides is unclear [5]. Understanding the constraints on the lithospheric thermal regime exerted by mantle dynamics has been considered the first step to study the high heat flow value (about 60 mWm −2 ) of the Internal Albanides [6].
A complex combination of various environmental factors controls the suitability of an area for the production of geothermal energy [7]. A geothermal resource is in fact part of a natural system in which geological features including rock type, diagenesis, mechanical characteristics of the rocks, faults, fractures, fluid chemistry, stress field and heat flow control key parameters such as the occurrence and spatial distribution of domains characterized by high porosity and high permeability (and related fluid circulation), vertical and lateral temperature gradients, and reservoir behavior during injection and production, which in turn is crucial for power plant efficiency. Therefore, a comprehensive picture of the geological setting and structural architecture of a potential geothermal site is fundamental for any site-specific, appropriate field development. For this reason, a prior geothermal suitability assessment is fundamental. This is commonly based on a series of exploration techniques often involving invasive inspections (e.g., well drilling), high costs and legal permissions [8]. Within this context, a reference model of the thermal structure at the crustal scale may represent a fundamental tool helping to save time and money during subsequent, local geothermal assessments. In this study, we compute surface heat flow, geotherms and isotherms along a balanced and restored regional geological section across northern Albania (Figure 1).
The resulting model provides a general framework of the thermal structure of this portion of the country, thus representing useful reference background for any further assessment of the geothermal energy systems of Albania. It is worth noting that, while prime geothermal systems occur in areas characterized by active volcanism/magmatism, the advancement of technology allow areas of moderate or even low heat flow to become of interest for geothermal exploration and exploitation.  The resulting model provides a general framework of the thermal structure of this portion of the country, thus representing useful reference background for any further assessment of the geothermal energy systems of Albania. It is worth noting that, while prime geothermal systems occur in areas characterized by active volcanism/magmatism, the advancement of technology allow areas of moderate or even low heat flow to become of interest for geothermal exploration and exploitation. Tectonic sketch map of the study area (modified from Roure et al. [10], Muceku et al. [11,12]). The black line A-A' shows the trace of the balanced and restored regional geological cross-section.
The study area ( Figure 2) is part of the outer Albanides fold and thrust belt [6,9,13] that derives from deformation of part of the former Adria continental margin, of African plate affinity [14,15]. This sector is characterized as NNW-SSE striking, ESE dipping thrusts subdividing the belt into a series of tectono-stratigraphic domains conventionally including, from east to west, the (i) Krasta, (ii) Kruja, (iii) Ionian and (iv) Sazani units [6]. The Sazani unit is continuous with the Apulia carbonate platform of the Apennines and Adriatic offshore; it is characterized by a late Triassic to Oligocene, shallow-water carbonate succession. This is overlain by a foreland basin succession mainly consisting of resedimented carbonates of Early Miocene-Pliocene age [16]. The Ionian basin started  [11,12]). The black line A-A' shows the trace of the balanced and restored regional geological cross-section.
The study area ( Figure 2) is part of the outer Albanides fold and thrust belt [6,9,13] that derives from deformation of part of the former Adria continental margin, of African plate affinity [14,15]. This sector is characterized as NNW-SSE striking, ESE dipping thrusts subdividing the belt into a series of tectono-stratigraphic domains conventionally including, from east to west, the (i) Krasta, (ii) Kruja, (iii) Ionian and (iv) Sazani units [6]. The Sazani unit is continuous with the Apulia carbonate platform of the Apennines and Adriatic offshore; it is characterized by a late Triassic to Oligocene, shallow-water carbonate succession. This is overlain by a foreland basin succession mainly consisting of resedimented carbonates of Early Miocene-Pliocene age [16]. The Ionian basin started Energies 2020, 13, 6028 4 of 19 to form in Triassic to Early Jurassic times by intracontinental rifting prior to the breakup that led eventually to seafloor spreading in the Tethys Ocean [14,17]. Further east, the Kruja carbonate platform and the Krasta unit face the so-called Internal (or inner) Albanides. The Krasta unit consists of Triassic siliciclastic strata and shallow-water limestones, Jurassic pelagic strata, Lower Cretaceous siliciclastic turbidites, Upper Cretaceous limestones and uppermost Cretaceous-Eocene siliciclastic turbidites [13]. The Internal Albanides consist of Triassic and Jurassic successions that experienced multiple deformation events associated with plate convergence [9].
The geological cross-section presented in this study is mainly showing Ionian basin units and stratigraphically overlying foreland basin deposits. The base of the Ionian basin succession consists of Triassic evaporites composed by intercalated gypsum, anhydrides and dolomitic limestone, resting on top of a crystalline basement [18]. Upward, the Lower Jurassic succession is formed by shallow water carbonates sedimented on the crest of tilted fault blocks overlapped by basin strata [19,20]. The Mesozoic to Paleogene succession is composed of pelagic limestones, marls and shales [19,20]. The overlying (up to 10 km thick) foreland basin succession fills the so-called Peri-Adriatic Depression. Turbiditic strata characterize the Oligo-Miocene part of the foreland basin succession [10,21]. Overlying Plio-Quaternary terrigenous sediments occur up to the Adriatic Sea bottom, also cropping out along the coastal areas of Albania [9,10] (Figure 2).

Active Tectonics
The study area is well known for a significant tectonic activity (with earthquakes M w ≤ 5.1) associated with active thrust faults, as also witnessed by historical seismicity [22][23][24] (Figure 3).
During 2019-2020, a seismic sequence affected the Albanian coastal area, north of Durrës, claiming 51 victims and causing severe damage. The INGV (Istituto Nazionale di Geofisica e Vulcanologia) recorded the main shock of the sequence, assigning to it a magnitude M w = 6.2. This event occurred on 26 November 2019, with an epicenter located 26 km northwest of Tirana, at a depth of 22 km [25].
From 21-09-2019 to 28-01-2020, 22 events with M w ≥ 3 occurred and were included in the INGV ISIDe catalogue [26]. They show a broad spatial distribution and a depth range of 10-36 km (mostly around 20 km). Available focal mechanisms show dominant thrust fault plane solutions, with mostly NNW-SSE striking nodal planes and sub-horizontal to gently dipping, mainly WSW-ENE trending P axes. These dates are coherent with the available seismic moment tensor solutions data by the RCMT catalogue, showing a dominant WSW-ENE trending maximum horizontal compression in the Durrës area and nearby Adriatic Sea. This seismic sequence provides important information not only on the seismotectonic setting of the region, but also on the crustal architecture of the outer Albanides and adjacent Adriatic crustal sector.
Geodetic calculations based on permanent GNSS Stations in the ETRF2000 "European fixed" reference frame [27] indicate that vector motion velocities in the western side of Albania are north-directed, while those in the central and eastern portions are south-directed, implying dextral strike-slip [28]. The projected motion vector velocities along the section trace A-A' in Figures 1 and 2 provide a south-westward velocity of 1.2 ± 0.7 mm yr −1 .  The epicenters are plotted using the location from the ISIDe [26], while magnitude (MW) and focal mechanism data are from the RCMT catalogue [25,29]. Historical earthquakes from 1000 to 1899 (blue squares) are selected by SHEEC catalogue [23,30] and from 1900 to 1984 by GFZ [22,31]. The co-seismic interferogram from Sentinel 1 [32] is shown as an overlay draped onto the digital elevation model [33]. . The epicenters are plotted using the location from the ISIDe [26], while magnitude (M W ) and focal mechanism data are from the RCMT catalogue [25,29]. Historical earthquakes from 1000 to 1899 (blue squares) are selected by SHEEC catalogue [23,30] and from 1900 to 1984 by GFZ [22,31]. The co-seismic interferogram from Sentinel 1 [32] is shown as an overlay draped onto the digital elevation model [33].

Materials and Methods
Aimed at modeling the thermal structure of the study area, a balanced and restored regional cross-section was built based on geological and geophysical datasets [6,10,12,34], as well as the new seismological dataset made available by the 2019-2020 seismic sequence [25]. The new geological cross-section provides pivotal information on the geometry and dimensional parameters of active thrusts in the frontal Albanides.
Using an analytical procedure, we calculated the geotherms for a series of pseudo-wells traced along the section. Interpolation of the computed geotherms allowed us to obtain relevant isotherms. We also defined the trend of the surface heat flow along the section by interpolating the computed surface heat flow for each pseudo-well.

Balanced and Restored Geological Cross-Section
Incorporating published geological and geophysical datasets, available seismic reflection data [6,10,12,34] and seismological information, we constructed a 120 km-long, regional ENE-WSW oriented cross-section, perpendicularly to the main structures of the area ( Figure 4).

Materials and Methods
Aimed at modeling the thermal structure of the study area, a balanced and restored regional cross-section was built based on geological and geophysical datasets [6,10,12,34], as well as the new seismological dataset made available by the 2019-2020 seismic sequence [25]. The new geological cross-section provides pivotal information on the geometry and dimensional parameters of active thrusts in the frontal Albanides.
Using an analytical procedure, we calculated the geotherms for a series of pseudo-wells traced along the section. Interpolation of the computed geotherms allowed us to obtain relevant isotherms. We also defined the trend of the surface heat flow along the section by interpolating the computed surface heat flow for each pseudo-well.

Balanced and Restored Geological Cross-Section
Incorporating published geological and geophysical datasets, available seismic reflection data [6,10,12,34] and seismological information, we constructed a 120 km-long, regional ENE-WSW oriented cross-section, perpendicularly to the main structures of the area ( Figure 4).  Figure 1 and Figure 2). The sedimentary cover structure is based on Frasheri et al. [6], Roure et al. [10] and Muceku et al. [12] (modified). The Moho is based on Grad and Tiire [34]. The section shows an earthquakes' hypocenter (some earthquakes clustering at a 10 km depth is an automatic location). No internal stratigraphy is shown for the Internal Albanides and Krasta zone. Pin lines used for cross-section restoration are shown.
The interpreted seismic profile by Roure et al. [10] provides constraints on the central and eastern portion of the geological section in which the seismic line shows the structures involving the Kruja carbonate platform and stratigraphically overlying foredeep clastic succession. On the eastern part of Figure 4, the Messinian unconformity truncates a W-transported thrust stack. In the footwall  Figures 1 and 2). The sedimentary cover structure is based on Frasheri et al. [6], Roure et al. [10] and Muceku et al. [12] (modified). The Moho is based on Grad and Tiire [34]. The section shows an earthquakes' hypocenter (some earthquakes clustering at a 10 km depth is an automatic location). No internal stratigraphy is shown for the Internal Albanides and Krasta zone. Pin lines used for cross-section restoration are shown.  Figure 4, the Messinian unconformity truncates a W-transported thrust stack. In the footwall of the latter, a low-displacement basement thrust is interpreted in the section [10]. A fish-tail geometry characterizes the upper and frontal portion of this structure. Two exploration wells penetrated two reverse faults-bounded pop-ups in the central part of the section. These structures are associated with Plio-Quaternary major faults extending at depth into the underlying Ionian basin strata. Roure et al. [10] interpreted the main thrust defining the eastern pop-up as a steep fault down to basement depths. The next pop-up is also likely to represent a shallow structure linked with a basement fault at depth. Frasheri et al. [6] showed an anticline further west, imaged by offshore geophysical datasets. This structure appears to represent a hanging-wall anticline related with the Albanides frontal thrust. This consists in a moderately dipping basement fault becoming gently dipping at the hypocentral depth of the M w = 6.2 main event of the 2019-2020 seismic sequence. This fault, representing a low-displacement thrust offsetting the Ionian basin succession, is interpreted to detach along a middle crustal decollement. Frasheri et al. [6] also show several extensional faults in the Sazani carbonate platform and associated paleoslope based on offshore seismic data. These inactive faults are associated with the extensional processes affecting the fore-bulge and the outer portions of the flexural foreland basins, and/or with the pre-orogenic inherited architectures of the Apulian platform paleo-slope.
The cross-section was line-length balanced using the Permian-Triassic level as a key marker for line-length measurements. Restoration was carried out using Move software (Petroleum Experts) for a segment in the footwall of the Kruja Thrust to a pin line placed in the foreland ( Figure 4). Fold restoration was carried out applying the flexural slip algorithm. Thrust restoration was carried out by using the fault parallel flow algorithm, most suitable for reverse fault restoration (also preserving bed area and length). The cumulative length of the Permian-Triassic level is shown in the restored cross-section of  of the latter, a low-displacement basement thrust is interpreted in the section [10]. A fish-tail geometry characterizes the upper and frontal portion of this structure. Two exploration wells penetrated two reverse faults-bounded pop-ups in the central part of the section. These structures are associated with Plio-Quaternary major faults extending at depth into the underlying Ionian basin strata. Roure et al. [10] interpreted the main thrust defining the eastern pop-up as a steep fault down to basement depths. The next pop-up is also likely to represent a shallow structure linked with a basement fault at depth. Frasheri et al. [6] showed an anticline further west, imaged by offshore geophysical datasets. This structure appears to represent a hanging-wall anticline related with the Albanides frontal thrust. This consists in a moderately dipping basement fault becoming gently dipping at the hypocentral depth of the Mw = 6.2 main event of the 2019-2020 seismic sequence. This fault, representing a low-displacement thrust offsetting the Ionian basin succession, is interpreted to detach along a middle crustal decollement. Frasheri et al. [6] also show several extensional faults in the Sazani carbonate platform and associated paleoslope based on offshore seismic data. These inactive faults are associated with the extensional processes affecting the fore-bulge and the outer portions of the flexural foreland basins, and/or with the pre-orogenic inherited architectures of the Apulian platform paleo-slope. The cross-section was line-length balanced using the Permian-Triassic level as a key marker for line-length measurements. Restoration was carried out using Move software (Petroleum Experts) for a segment in the footwall of the Kruja Thrust to a pin line placed in the foreland ( Figure 4). Fold restoration was carried out applying the flexural slip algorithm. Thrust restoration was carried out by using the fault parallel flow algorithm, most suitable for reverse fault restoration (also preserving bed area and length). The cumulative length of the Permian-Triassic level is shown in the restored cross-section of Figure 5

Analytical Procedure for Thermal Modeling
The thermal model is produced on the basis of the study of Molnar et al. [35], modified after Candela et al. [36]. Three sources of heat affect the thermal structure and the resulting heat flow measured at the surface ( ): (i) the heat produced by the mantle, considered here as flowing

Analytical Procedure for Thermal Modeling
The thermal model is produced on the basis of the study of Molnar et al. [35], modified after Candela et al. [36]. Three sources of heat affect the thermal structure and the resulting heat flow measured at the surface (Q S ): (i) the heat produced by the mantle, considered here as flowing upward from the Moho discontinuity (Q m ), (ii) the heat production rate (H) due to radioactive decay in the crust, and (iii) the frictional heating produced by thrust faulting. The analytical procedure was applied in a series of pseudo-wells traced on the balanced geological cross-section ( Figure 6) taking into account present-day topography (Figure 1). For each pseudo-well we calculated the geotherm. The results were interpolated through a polynomial equation to obtain an overall picture of the surface heat flow density and isotherms from the offshore to eastern Albania.
Energies 2020, 13, x FOR PEER REVIEW 8 of 21 upward from the Moho discontinuity ( ), (ii) the heat production rate ( ) due to radioactive decay in the crust, and (iii) the frictional heating produced by thrust faulting. The analytical procedure was applied in a series of pseudo-wells traced on the balanced geological cross-section ( Figure 6) taking into account present-day topography (Figure 1). For each pseudo-well we calculated the geotherm. The results were interpolated through a polynomial equation to obtain an overall picture of the surface heat flow density and isotherms from the offshore to eastern Albania. In the analysis, we refer to the following parameters: depth ( ), thickness of the sedimentary cover (ℎ ), thickness of the basement (ℎ ), temperature due to mantle heat flow ( ), radiogenic heat ( ) and frictional heat ( ), respectively. The geotherms were obtained separately for the sedimentary cover (for a depth ≤ ℎ ) and for the crystalline upper crust (ℎ < ≤ ℎ, with ℎ = ℎ + ℎ ). Concerning the latter, frictional heating ( ) was taken into account in addition to the two main heat sources ( and ) associated with the new, perturbed conductive state. For the upper layer, the increment of radiogenic heat produced by basement thrusting at depth was computed at the bottom of the sedimentary cover, adding to the previous equilibrium status ( ). Frictional heating ( ) in the sedimentary cover was also taken into account. The adopted analytical approach is described by the following equations: where is the temperature at the surface and (ℎ , ) is the temperature at the bottom of the sedimentary cover. Each time-dependent term in these equations has a generic form of the type: where ( ) is the final temperature at the new state of equilibrium ( → ∞) and is the thermal diffusivity coefficient ( 1.0 mm 2 s −1 for all time-dependent terms). All factors for each specific expression are shown in Table 1, while the used parameters are listed in Table 2. In the analysis, we refer to the following parameters: depth (z), thickness of the sedimentary cover (h C ), thickness of the basement (h B ), temperature due to mantle heat flow (T MH ), radiogenic heat (T RH ) and frictional heat (T FH ), respectively. The geotherms were obtained separately for the sedimentary cover (for a depth z ≤ h C ) and for the crystalline upper crust (h C < z ≤ h, with h = h C + h B ). Concerning the latter, frictional heating (T FH ) was taken into account in addition to the two main heat sources (T MH and T RH ) associated with the new, perturbed conductive state. For the upper layer, the increment of radiogenic heat produced by basement thrusting at depth was computed at the bottom of the sedimentary cover, adding to the previous equilibrium status (T DH ). Frictional heating (T FH ) in the sedimentary cover was also taken into account. The adopted analytical approach is described by the following equations: where T S is the temperature at the surface and T DH (h C , t) is the temperature at the bottom of the sedimentary cover. Each time-dependent term in these equations has a generic form of the type: A n sin(a n z) −a 2 n Kt , X ≡ D, M, R, F where T f (z) is the final temperature at the new state of equilibrium ( t → ∞ ) and K is the thermal diffusivity coefficient (1.0 mm 2 s −1 for all time-dependent terms). All factors for each specific expression are shown in Table 1, while the used parameters are listed in Table 2. Table 1. List of terms included in the equations adopted in the analytical study.
A value of Q m = 10 mWm −2 was assumed at the lower boundary of the model, according to the Moho heat flow trend obtained for the Albania offshore by Cermak [40]. We used the values obtained by Dragoni et al. [41] for the basement of the Apulian foreland crust for the heat related with radioactive decay in the crust (H B = 3.15 mWm −3 ) and for the distance of its decrease with depth (H = 8 km). Instead, we assumed different values for the radioactive decay-related heat (H B < H B ), i.e., 1.6 mWm −1 for the internal part of the Internal Albanides crystalline basement (pseudo-wells number 4 and 5 in Figure 6) [5,11], considering the observation by Frasheri and Frasheri [42] that "the granites of the crystalline basement, with the radiogenic heat generation, represent the heat source". This takes into account the very low radiogenic heat produced by ophiolitic rocks, as well as the increasing surface heat flow measured in the ophiolitic belt being mainly produced by heat rising from depth [42].
In the computation of the surface heat flow, in addition to contributions of Q m and the radiogenic heat, we also added the frictional heat due to overthrusting: Figure 4 shows our 120 km-long, balanced and restored cross-section. Restoration of the regional balanced section ( Figure 5) allowed us to calculate the cumulative shortening, resulting in 9.10 km. Fault 1, which nucleated the main shock of the 2019-2020 seismic sequence, is responsible for~13% (1.20 km) of the cumulative displacement; fault 2 for~8% (0.75 km); fault 3 for~9% (0.80 km), fault 4 for~19% (1.69 km), and fault 5 for~36% (3.29 km). In total, the shortening produced by faulting is 7.73 km, i.e.,~85% of the total. The remaining amount of shortening is due to folding, mainly associated with the south-westernmost frontal anticline (Figure 4). Our balanced and restored regional geological section highlights the basement-involved thrusting style characterizing the outer portion of the Albanides fold and thrust belt at the latitude of Durrës (~41 20 N). The cross-section (Figure 4) was based on previous published sections, geophysical dataset datasets [6,10,12,34] and seismological information [25,26]. Various studies, e.g., [14], documented a different structural style more to the south, where thin-skinned thrusting of the sedimentary cover was promoted by the occurrence of an efficient decollement level represented by the Triassic evaporites (which also form large diapirs). In the study sector, instead, the Triassic units appear to have a negligible effect and there is no basement-cover decoupling, most probably due to their reduced thickness and decreasing amount of evaporite. The change in tectonic style from thin-skinned thrusting involving the sedimentary cover alone to basement-involved thrusting seems to be controlled by NE-SW striking transfer zones (e.g., Shkoder-Peje, Vlora-Elbasan) interpreted to have reworked pre-existing basement structures [43,44]. Towards the foreland, the inherited continental margin architecture most probably also controlled the large accommodation space for the 10 km thick foreland basin succession hosted in the Peri-Adriatic Depression, as documented in other fold and thrust belts characterized by marked along-strike variations of foreland basin geometry, accommodation space and related sediment accumulation [45].

Balanced Cross-Section
The restored section ( Figure 5) shows how the cumulative shortening produced by basement faults is low (7.730 km), although it represents~85% of the total cumulative shortening (the remaining shortening being associated with folding). According to the seismic evidence provided by Roure et al. [10], the shortening related to basement-involved thrusting was accumulated in the last 5 Ma. The resulting shortening rates are consistent with recent geodetic data of permanent GNSS Stations in the ETRF2000 "European fixed" reference frame [36], which yield a southward motion of 1.2 ± 0.7 mm yr −1 projected into our section of Figure 4.

Thermal Modeling
The analytical procedure described in Section 3.2 allowed us to calculate a geotherm (Figure 7) for each pseudo-well located in the balanced cross-section ( Figure 6), fixing T S = 15 • C and considering also the topography for pseudo-wells 4, 5 and 6. Consequently, the altitude causes a different value of temperature at the level of our fixed reference level (z = 0) with respect to T S = 15 • C; e.g., in the case of pseudo-well number 6, where the surface elevation is about 2.7 km, the temperature to z = 0 is ca.
The presence of the thrust produces an additional heat source, which causes an increase in the thermal gradient above to the depth of the thrust fault. Instead, at greater depths (i.e., below the thrust) the thermal gradient is less because of the absence of the frictional heat. Therefore, all the geotherms show a temperature increase above the depth at which the presence of the thrust produces a change of the curve inclination. Moreover, the thermal gradient of geotherms 1, 2, 3 (Peri-Adriatic Depression) and 6 ranges from 14.6 to 16.8 mKm −1 , these values being similar to those obtained by Frasheri et al. [6] for the Pre-Adriatic Depression (from 15 to 21.3 mKm −1 ). Instead, geotherms 4 and 5 show a different temperature gradient (9.7-10.5 mKm −1 ) produced by the minor radiogenic heat (H B ) imposed, following Muceku et al. [12], taking into account Cermak et al. [5].
Megna et al. [46] checked the sensitivity of this type of analytical procedure calculating that a change of about 10% of slip rate or thrust depth would produce a variation of frictional heat (T FH ) able to modify the resulting temperature by a value always below 2%. This type of analytical procedure was used also in subsequent works by Basilici et al. [47,48], confirming a similar limited error. The computed surface heat flux Q s shown in Figure 8a (orange line) was interpolated using a fifth-order polynomial equation, which represents the best fitting analytical curve for the geothermal values. In the western area, relative high values are located between pseudo-wells 2 and 3 due to the frictional heat produced by faults. The lower values are reached between pseudo-wells 4 and 5 due to a lower radiogenic heat (H B ) produced by the basement. The highest values of Q s are reached to the east due to the radiogenic heat (H B ), the presence of a high number of faults, and the highly conductive and slightly radioactive ophiolite of the inner Albanides.
The isotherms shown in Figure 8b were interpolated using a third-order polynomial equation, also in this case representing the best fit for calculated values. They provide an overall NE-ward increasing temperature pattern, however, with a concavity defining minimum values in correspondence of pseudo-well 3. In the Adriatic sector, isotherms rise up and become closer to each other. The presence of the thrust produces an additional heat source, which causes an increase in the thermal gradient above to the depth of the thrust fault. Instead, at greater depths (i.e., below the thrust) the thermal gradient is less because of the absence of the frictional heat. Therefore, all the geotherms show a temperature increase above the depth at which the presence of the thrust produces a change of the curve inclination. Moreover, the thermal gradient of geotherms 1, 2, 3  conductive and slightly radioactive ophiolite of the inner Albanides.
The isotherms shown in Figure 8b were interpolated using a third-order polynomial equation, also in this case representing the best fit for calculated values. They provide an overall NE-ward increasing temperature pattern, however, with a concavity defining minimum values in correspondence of pseudo-well 3. In the Adriatic sector, isotherms rise up and become closer to each other. The geotherms (Figure 7) obtained for a series of six pseudo-wells traced on the section ( Figure  6) effectively display the thermal structure of the northern sector of outer Albanides and adjacent Adriatic sector, including isotherms and surface heat flow (Figure 8). The south-westernmost portion of the thermal model shows how the isotherms above the main frontal thrust (i.e., those with temperature ≤ 250 °C) are closely spaced and roughly equidistant, whereas the 300 °C isotherm lying at a depth greater than the thrust is more distant and significantly depressed due to the different thermal gradient below the thrust (Figure 7). This feature confirms that frictional heating produced by faulting strongly influences the thermal structure of fold and thrust belts, raising the isotherms [47,48]. The geotherms (Figure 7) obtained for a series of six pseudo-wells traced on the section ( Figure 6) effectively display the thermal structure of the northern sector of outer Albanides and adjacent Adriatic sector, including isotherms and surface heat flow ( Figure 8). The south-westernmost portion of the thermal model shows how the isotherms above the main frontal thrust (i.e., those with temperature ≤ 250 • C) are closely spaced and roughly equidistant, whereas the 300 • C isotherm lying at a depth greater than the thrust is more distant and significantly depressed due to the different thermal gradient below the thrust (Figure 7). This feature confirms that frictional heating produced by faulting strongly influences the thermal structure of fold and thrust belts, raising the isotherms [47,48].
To improve the description of our results, the model was compared with previous studies on the thermal structure [6] along a section located close to our geological cross-section ( Figure 9). Heat flow curves show a similar trend, with lower values within the Kruja area (85-100 km along section) and higher values on the Peri-Adriatic Depession and Internal Albanides. The computed data show overall higher values with respect to those obtained by Frasheri et al. [6], the comparison yielding a root main square (rms) value of 10.3 mWm −2 . The temperature isolines trend is generally in agreement with that presented by Frasheri et al. [6] (Figure 9b). Some differences (in term of values) occur because we considered a different Moho depth [34] and the influence of topography, which is negligible in the Adriatic sector but becomes significant eastward, as visible by the starting point of the geotherms in Figure 7. Nevertheless, there is a good fit (rms = 1.6 • C) between our geotherm calculated for the Peri-Adriatic Depression (pseudo-well 2; Figure 10) and the temperature profile provided by Frasheri [49] for a borehole located ca. 100 km along the strike to the south. root main square (rms) value of 10.3 mWm −2 . The temperature isolines trend is generally in agreement with that presented by Frasheri et al. [6] (Figure 9b). Some differences (in term of values) occur because we considered a different Moho depth [34] and the influence of topography, which is negligible in the Adriatic sector but becomes significant eastward, as visible by the starting point of the geotherms in Figure 7. Nevertheless, there is a good fit (rms = 1.6 °C) between our geotherm calculated for the Peri-Adriatic Depression (pseudo-well 2; Figure 10) and the temperature profile provided by Frasheri [49] for a borehole located ca. 100 km along the strike to the south.   [6]. (b) Thermal structure of the outer Albanides and adjacent Adriatic sector from this study, compared with isotherms by Frasheri et al. [6].
Energies 2020, 13, x FOR PEER REVIEW 16 of 21 Figure 10. Temperature profile for the Peri-Adriatic Depression from Frasheri [49] modified and compared with the geotherm obtained in this study for a similar structural setting along strike (pseudo-well 2, located in Figure 6). The top of the carbonates (green) is shallower with respect to that shown in the geological section of this study (Figure 4) due to the northward plunge of the structures.
values obtained in this study display higher values with respect to those of Frasheri et al. [6], however showing a similar shape of the curve (Figure 9a). The basement represents the main radiogenic heat source; its composition may vary strongly [12]. For this reason, we suggest that minimum values of located in the Tirana area [5,6,42,49], corresponding with the NE part of the section, are due to a lower radiogenic heat ( ′ ) produced by the basement in that zone. Figure 10. Temperature profile for the Peri-Adriatic Depression from Frasheri [49] modified and compared with the geotherm obtained in this study for a similar structural setting along strike (pseudo-well 2, located in Figure 6). The top of the carbonates (green) is shallower with respect to that shown in the geological section of this study (Figure 4) due to the northward plunge of the structures.

Potential Exploitation of the Geothermal Energy
Q s values obtained in this study display higher values with respect to those of Frasheri et al. [6], however showing a similar shape of the curve (Figure 9a). The basement represents the main radiogenic heat source; its composition may vary strongly [12]. For this reason, we suggest that minimum values of Q s located in the Tirana area [5,6,42,49], corresponding with the NE part of the section, are due to a lower radiogenic heat (H B ) produced by the basement in that zone.

Potential Exploitation of the Geothermal Energy
Minimum values of Q s in the surroundings of Tirana and the reconstructed 2D thermal structure would suggest less favorable conditions for exploitation of geothermal energy, besides the direct use (Borehole Heat Exchanger-Geothermal Heat Pump systems). Nevertheless, the "Kruja geothermal zone" [42], extending for 180 km with several hot spring manifestations, has an NNW-SSE trend and it is partially located NW and W of Tirana. This emphasizes the importance of tectonic features in driving hot fluids to the surface with respect to the general thermal structure. In fact, the trend of the "Kruja geothermal zone" corresponds to a major imbricate thrust fan involving Mesozoic-Paleogene limestones, partially sealed by Oligo-Miocene turbidites and by Plio-Quaternary terrigenous sediments overlying a Messinian unconformity (Figure 4). The role of thrust zones in fluid migration at the orogen scale is well known [50][51][52][53]. It may be envisaged that meteoric water infiltrating within outcropping fractured units of the thrust belt reaches depths of 7-8 km, where it reaches maximum temperatures >100 • C. The mineralized hot water is then channeled along the thrust zones of the Kruja imbricate fan ( Figure 11). Where the hot fluids can reach the surface, they give rise to a belt of mineralized thermal water springs [49]. On the other hand, where the carbonate anticlinal traps are sealed by Oligocene to Neogene-Quaternary terrigenous units, a drillable geothermal play may occur.
Energies 2020, 13, x FOR PEER REVIEW 17 of 21 fact, the trend of the "Kruja geothermal zone" corresponds to a major imbricate thrust fan involving Mesozoic-Paleogene limestones, partially sealed by Oligo-Miocene turbidites and by Plio-Quaternary terrigenous sediments overlying a Messinian unconformity (Figure 4). The role of thrust zones in fluid migration at the orogen scale is well known [50][51][52][53]. It may be envisaged that meteoric water infiltrating within outcropping fractured units of the thrust belt reaches depths of 7-8 km, where it reaches maximum temperatures >100 °C. The mineralized hot water is then channeled along the thrust zones of the Kruja imbricate fan (Figure 11). Where the hot fluids can reach the surface, they give rise to a belt of mineralized thermal water springs [49]. On the other hand, where the carbonate anticlinal traps are sealed by Oligocene to Neogene-Quaternary terrigenous units, a drillable geothermal play may occur. The higher surface heat flow characterizing the Internal Albanides, associated with more intense radiogenic heat generation related with basement composition, can also suggest a favorable condition for potential geothermal development, particularly in view of EGS. Thanks to modern technology, low enthalpy geothermal energy can be exploited in populated areas with good success and better yields especially in coastal areas. In addition, the areas with hot thermal spring manifestations are suitable for balneotherapy, greenhouses, and aquaculture/aquatic farming purposes. In the vicinity of inhabited centers, the possibility of developing geothermal binary systems for electricity production, together with heating and cooling and integrated with other renewable energy sources, is not excluded if the fluid temperature conditions (90°-100°) at depths not greater than 1000-1500 m can be verified with dedicated drillings. Figure 11. Detail of the regional cross-section, showing a conceptual model of fluid migration pathways controlling the "Kruja geothermal zone". The higher surface heat flow characterizing the Internal Albanides, associated with more intense radiogenic heat generation related with basement composition, can also suggest a favorable condition for potential geothermal development, particularly in view of EGS. Thanks to modern technology, low enthalpy geothermal energy can be exploited in populated areas with good success and better yields especially in coastal areas. In addition, the areas with hot thermal spring manifestations are suitable for balneotherapy, greenhouses, and aquaculture/aquatic farming purposes. In the vicinity of inhabited centers, the possibility of developing geothermal binary systems for electricity production, together with heating and cooling and integrated with other renewable energy sources, is not excluded if the fluid temperature conditions (90 • -100 • ) at depths not greater than 1000-1500 m can be verified with dedicated drillings.

Conclusions
The results of this study provide the first comprehensive account of the thermal structure of the northern sector of the outer Albanides and adjacent Adriatic sector obtained by an analytical methodology. The thermal model is portrayed along a balanced and restored regional geological cross-section produced using published and unpublished subsurface datasets, integrated with new seismological data related with the 2019-2020 seismic sequence. The thermal model of the study area, obtained using an analytical procedure, was compared with previous studies to improve the description of this geologically important area. We highlighted a zone (Tirana area) in which the main radiogenic heat source, represented by the basement, produces a minor amount of radiogenic heat, probably due to its varying composition [12]. The occurrence of the "Kruja geothermal zone", partially overlapping this latter area, points to the fundamental role exerted by the structural setting in the rise of hot fluids toward the surface. On the other hand, the inner Albanides to the E and the Adriatic sector to the W display a thermal setting strongly influenced by frictional heating produced by faults, similarly to what was reported for other active collisional orogens [47,48].
The thermal model obtained in this study may be used as background information on the geothermal suitability of different regions of northern Albania, prior to the potential application of expensive exploration techniques and well drilling. Our results may contribute to lowering mining risks, thereby allowing geothermal companies to save investment time and money. Furthermore, they can help policy makers in defining more effective energy management strategies. The publication of thermal models such as the present one may also aid to inform inhabitants about the geothermal energy potential of the areas they reside in, thus allowing the citizenship to benefit from enhanced safety and transparency of future decisions, and geothermal companies to make more secure investments.
Taking into account the thermal structure of the crustal sector analyzed in this study, further development of the geothermal resource to produce renewable energy in Albania is advisable. In particular, there seem to be suitable conditions for a wider dissemination of low-enthalpy small plants for heating/cooling of public and private buildings (direct use of heat). In specific areas, dedicated studies should be aimed at verifying the possibility of binary systems plants (medium enthalpy), perhaps integrated with other renewable energies (solar, biomass, wind).