Central Anatolia is a region of lithospheric thinning (e.g., Piromallo and Morelli, 2003; Fichtner et al., 2013) hosted between the Western Anatolia Extentional Province and the Eastern Anatolia Shortening Province (e.g., Göğüş et al., 2017, and references therein) (Figure 1). The tectonics of Central Anatolia are dominated by two active NW-SE trending strike-slip fault systems: 1) the left-lateral Ecemiş fault zone and 2) the transtensive dextral Tuz Gölü fault zone (TGFZ) (Toprak and Göncöoḡlu, 1993). In addition, sets of older normal faults with a NE-SW trend have been identified buried beneath Pleistocene volcanic products (Toprak and Göncöoḡlu, 1993; Kürçer and Gökten, 2014). The 200 km long (Kürçer and Gökten, 2014) TGFZ marks the western border of the Central Anatolian Volcanic Province (CAVP). Seismicity along the fault is common, with an average earthquake focal depth of 10 km, and a moment magnitude M _w ≤ 5.2 (Kürçer and Gökten, 2014). Block modelling across Central Anatolia by Aktuğ et al. (2013) indicates a normal component of the slip rate along the TGFZ of about 1 mm yr⁻¹, in addition to a right-lateral slip rate of about 5 mm yr⁻¹ along the middle and/or southern part of the TGFZ (Aktuğ et al., 2013).

The TGFZ is composed of fault segments with lengths of between 4 and 33 km (Kürçer and Gökten, 2012). With a length of 27 km, the NW-SE trending Akhisar-Kiliç fault segment (AKFS) is one of the most important segments of the TGFZ due to its length and morphotectonic features which indicate a total vertical displacement of 268 m since the lower Pliocene (Kürçer and Gökten, 2012).

Hasan Dağ is a polygenetic, twin peaked stratovolcano in the CAVP (Kuscu and Geneli, 2010) composed of Big Mount Hasan (BMH) at 3252 m a.s.l. and Small Mount Hasan (SMH) at 3069 m a.s.l. (Figure 1). Edifice formation occurred on the site of the extinct Keçikalesi volcano by three, distinctive evolutionary phases (Paleo-, Meso-, Neo-Hasan Dağ; < 7.21 ± 0.9 Ma, 1-0.1 Ma, and < 0.1 Ma, respectively) (Aydar, 1992; Aydar and Gourgaud, 1998). All evolutionary phases have temporal patterns synonymous with ignimbrite volcanism: dome lava extrusion preceding extensive volcanoclastic deposits and structural failure of the vent leading to caldera formation (Deniel et al., 1998). Magmas at Hasan Dağ have evolved from tholeiitic to calc-alkaline but are consistently dominated by andesites and rhyolites with subordinate basalts (Deniel et al., 1998).

Hasan Dağ’s edifice is composed of pyroclastic deposits, lava flows and lava domes (Friedrichs et al., 2020), with the most recent pumice-forming eruption dated at 8.98 ± 0.64 ka (Schmitt et al., 2014) and a lava dome dated at < 6 ka (Aydar and Gourgaud, 1998). A long-term eruptive flux for Hasan Dağ is estimated at between 0.1 and 0.3 km³ ka⁻¹ (Friedrichs et al., 2020). Elevated levels of H₂S and CO₂ gas at the volcano’s summit and fumarolic degassing on its western flank with temperatures of up to 68.7°C and CO₂ levels ∼10,000 ppm are indicative of the presence of an active magmatic system (Dİker et al., 2018; Ulusoy et al., 2020; Aydar et al., 2021).

Mid-crustal low-seismic velocity layers beneath the volcanic province may represent magma reservoirs feeding active surface volcanism (Abgarmi et al., 2017). The electrical conductivity structure beneath Hasan Dağ highlights a conductive zone centred at 5 km depth beneath the plateau and is interpreted to represent the volcano’s upper-crustal magma reservoir (Tank and Karaş, 2020).

The fault segment (AKFS) transects the north-eastern flank of the Hasan Dağ volcanic complex (Kürçer and Gökten, 2014) and intersects the sinistral, NE-SW trending Karaman-Niğde shear zone (KNSZ) (Aydar and Gourgaud, 1998; Şengör and Zabcı, 2019) (Figure 1). Quaternary volcanics of Hasan Dağ have been cut by faulting along the AKFS (Kürçer and Gökten, 2014). The segment has an average Quaternary vertical slip rate of 0.05 mm yr⁻¹ (Kürçer and Gökten, 2012). Dextral offsets of Hasan Dağ lavas permit deducing average right-lateral slip rates of 0.58–3.35 mm yr⁻¹ over the past 2 Ma (Krystopowicz et al., 2020). Paleoseismic and geochronological evidence (Kürçer and Gökten, 2012) of two Holocene earthquakes on the AKFS at 10 ± 0.41 ka BP and 5.450 ka±0.16 BP provide a recurrence interval of 4.55 ± 0.45 ka. Elevated footwall uplift rates towards the south-eastern end of the TGFZ may pertain to interactions between the fault zone and the Hasan Dağ magmatic system (Krystopowicz et al., 2020). Seismicity in proximity to Hasan Dağ includes hypocentral depths that are deeper ( > 10 km) than common for the region (Aydar et al., 2021), and may perhaps be indicative of a volcanogenic origin (Kürçer and Gökten, 2014).

To investigate volcano-tectonic linkages between the Tuz Gölü fault segment and Hasan Dağ we numerically solve for stresses imparted by slip on the fault segment (the source) to the surrounding crust and faults/fractures (the receivers) using Finite Element Analysis (COMSOL Multiphysics v5.4). The 3D model (Figure 2) is set up such that horizontal distances were −45 km ≤ x ≤ 45 km and −75 km ≤ y ≤ 75 km, while the model extends to a depth z of 36 km reflecting the depth of the Moho beneath the plateau surrounding Hasan Dağ (Tezel et al., 2013). The plateau sits at 1000 m above sea level with Hasan Dağ’s edifice protruding above to a maximum prominence of 2253 m. Geospatial data of the region were obtained from a Digital Elevation Model (Aktuğ et al., 2013) based on the 3 arc-seconds NASA Shuttle Radar Topographic Mission (SRTM) data (Version 2; https://www2.jpl.nasa.gov/srtm/).

Model boundary condition are: 1) a free surface boundary constraint applied to the upper surface of the model, 2) a fixed constrained applied to the lower boundary, and 3) roller constraints at the lateral surfaces. The TGFZ and KNSZ are implemented as a set of identity pairs with continuity constraints applied and have either slip velocities or finite displacements prescribed depending on the nature of the study. The KNSZ has a no-slip boundary condition attributed to it. Both fault zones were modelled as straight, continuous features running through the entire modelling domain. While the TGFZ is modelled as a subvertical right-lateral fault trending 34°NW and dipping 68°SW from the horizontal, the KNSZ is modelled as vertical and trending at 55°NE (Kürçer and Gökten, 2014).

Using the upper-crustal high electrical conductivity anomaly beneath the volcano as a proxy (Tank and Karaş, 2020), we implement the volcano’s magma reservoir by an oblate ellipsoid with a semi-minor axis of 1 km and semi-major axes of 1.8 km. The modelled reservoir contains six (a–f) point probes (Figure 2). a is at the upper most point on the reservoir while b is at the lower-most point. Both are 1 km from the reservoir’s centre. Points a and b have tangents that are in the xy plane. c, d, e and f are on the side of the magma reservoir at greatest distance from the centre (1.8 km). Points c and d have tangents that are parallel to the fault plane of the AKFS while e and f have tangents that are perpendicular to it. The magma reservoir is located towards the south-eastern end of the fault plane.

A typical model contains a structured mesh with ∼ 2.5 ×10⁵ elements solving for ∼10⁶ degrees of freedom. The mesh coarsens with distance from the fault plane, magma reservoir, and edifice.

Surface rupture lengths of 5, 10, 15, 20, 25, and 27 km were used to parameterise the partial or complete failure of the AKFS. To derive the moment magnitude M _w of the consequential earthquakes, we use the empirical relationship defined by Wells and Coppersmith (1994): M w = a + b l o g S R L (5) where SRL is the surface rupture length in km. a and b are 5.16±0.13 and 1.12±0.08, respectively, for strike-slip faulting, and 4.86±0.34 and 1.32±0.26, respectively for normal slip.

The fault rupture area (RA) is calculated from: l o g R A = a + b M w (6)

If the faulting regime is strike-slip then a and b are −3.42±0.18 and 0.9±0.03, respectively, while for normal slip they are −2.87±0.50 and 0.82±0.08, respectively (Wells and Coppersmith, 1994).

Fault rupture widths (RW; in km) in down-dip direction for given M _w were derived from: l o g R W = a + b M w (7) where a and b are −0.76±0.12 and 0.27±0.02, respectively, for strike-slip faulting, and 1.14±0.28 and 0.35±0.05, respectively for normal slip (Wells and Coppersmith, 1994).

See Table 5 for the full parameterisation.

We benchmarked the numerical solutions of Coulomb failure stress changes against analytical solutions from the Coulomb 3.4 code (Toda et al., 2011) for the case of an elastic and mechanically homogeneous medium. Coulomb 3.4 outputs ΔCFS in 2D space, while our models solve in 3D. The simplest way to benchmark the 2D and 3D solutions is to solve for ΔCFS in the xy, xz and yz planes in our models and to identify the plane of both maximum and minimum principal stresses and hence of maximum shear stresses. Due to the dominance of right-lateral displacement along the sub-vertical AKFS, the Coulomb failure stress changes calculated in the xy plane in our 3D models provided the best fit to the Coulomb 3.4 stress solutions. As a result, all ΔCFS solutions from our numerical models are reported in the xy plane such that ΔCFS = ΔCFF _xy . See Supplementary Material for further details on model benchmarking.

The mechanical structure of the crust is parameterised from seismic p-wave velocity data. Data from Salah et al. (2014) are used for depths > 3 km, while in the absence of reliable shallow ( < 3 km depth) crustal data for the CAVP and Hasan Dağ, the upper 3 km of the crust are parameterised using data from the Cascades. We chose seismic velocity data from Mt. St. Helens Volcano Kiser et al. (2016) due to close similarities in eruption style, eruptive products and edifice formation at both volcanoes and very similar upper-crustal seismic velocities (e.g., at 3 km depth V _p and V _p /V _s data in Central Anatolia are only 1.5% larger than in the Cascades). Note that the derived Coulomb Failure stress changes reported in this study for depth > 4 km are insensitive to the implemented mechanical structure < 3 km depth. The crustal density (ρ) distribution was obtained from p-wave velocities via the following polynomial regression fit (Brocher, 2005): ρ = 1.6 612 V p − 0.4 721 V p 2 + 0.0 671 V p 3 − 0.0 043 V p 4 + 0.000 106 V p 5 (8)

Poisson’s ratio (ν) was calculated as a function of the p- and s-wave velocity ratios: ν = V p / V s 2 − 2 / V p / V s 2 − 1 (9)

Finally, the dynamic Young’s modulus (E) was calculated from: E d = V p 2 ν 1 + ν 1 − 2 ν / 1 − ν (10)

An empirical correction is applied to convert the dynamic modulus E _d to static modulus E _s . Following (Gudmundsson, 1988) the ratio of E _s /E _d used in this study is 0.5. Elastic properties attributed to Hasan Dağ’s edifice were taken from Montserrat (Young and Gottsmann, 2015) due to a lack of local data and broadly similar eruptive products of andesite volcanism at both volcanoes. See Table 3 for mechanical properties assigned to all model domains.

The moment magnitude (M _w ) is dimensionless, but can be linked to the seismic moment (M _o ) of an earthquake (Hanks and Kanamori, 1979): M w = 2 / 3 ⁡ log M o − 10.7 (11)

The above equation gives seismic moment in 10^–7 Nm. The seismic moment is a function of the stress drop Δσ experienced by the fault during an earthquake (Kanamori and Anderson, 1975): Δ σ = 2 M o / π R W 2 S R L (12)

Here we assume that the stress drop equates to the magnitude of stress required for the fault to slip; i.e., rupture occurs at ΔCFS ≥ Δ σ, and is further assumed to occur in the middle of the fault segment propagating outwards. We thus model the fault segment as transcurrent.

The cumulative displacement along a fault is a function of time (t), thus a time-dependent study is added to the model, whereby t = 0 is immediately after the most recent earthquake when the stress on the AKFS was at its minimum. This enables the interrogation of the temporal evolution of stress accumulation along fault segments with lengths equalling hypothetical future surface rupture lengths with a view to determine the timing of failure for given fault kinematics.

Changes in Coulomb failure stresses are modeled by 1) slip of the AKFS since its most recent earthquake (termed gradual slip hereafter) and 2) by a hypothetical earthquake caused by rupture of the AKFS. Quaternary (gradual) slip rates of 0.05 mm yr⁻¹ normal slip, and 0.58–3.35 mm yr⁻¹ right-lateral slip are applied to the AKFS until the calculated stress drop magnitude Δσ is reached. Instantaneous displacements vectors were then applied to the fault segment to model the resultant Coulomb failure stress change on receiver faults from the hypothetical earthquake. See also Supplementary Material for the parameterisation of Coulomb failure stress changes.

The average fault displacement AD due to each earthquake was calculated via l o g A D = a + b M w (13) where M in the moment magnitude of each slip event, a and b are constants equalling −6.32±0.61 and 0.9±0.09, respectively (Wells and Coppersmith, 1994). Slip as a result of the hypothetical earthquakes follows the rake of the gradual slip. Earthquakes following normal slip displacement of 0.053 6 mm yr⁻¹ and right-lateral slip of 0.58 mm yr⁻¹ events have a rake of −175°. Earthquakes following normal slip displacement of 0.053 6 mm yr⁻¹ and right-lateral slip of 3.35 mm yr⁻¹ events have a rake of −179°. The derived average displacements and slip components are given in Table 4.

The Coulomb failure stress change on Hasan Dağ’s magma reservoir is calculated for three different modes of crustal clamping and unclamping: 1) unclamping or clamping of magma pathways above the reservoir, 2) clamping or unclamping of pathways below the reservoir, 3) direct clamping of the reservoir walls. Six point probes on the magma reservoir are chosen to derive the stress change (Figure 2B). At these points, Coulomb failure stress change is quantified on the tangential plane to the reservoir or fracture surface (see Supplementary Material). Scenario 1 is calculated at point a both parallel and perpendicular to the strike of fault plane. Scenario 2 is calculated at point b also parallel and perpendicular to the strike of the fault. In both scenarios fractures are assumed vertical, dipping 90° from the horizontal. Scenario 3 is calculated at points c, d, e, and f. The temporal evolution of stress changes on the reservoir is calculated from the most recent earthquake up to just before each hypothetical earthquake. The resultant stress change due to each earthquake was quantified by the prescribed slip on the fault plane.

Table 5 shows the moment magnitudes (M _w ) of the hypothetical earthquakes along the AKFS, and their associated dimensions for either strike-slip or normal slip kinematics. The fault kinematics have little effect on the moment magnitudes for equivalent surface rupture lengths (SRLs). Moment magnitudes of strike-slip faulting match those by normal slip faulting for SRLs ≥10 km. The difference in M _w for SRL = 5 km is <3%. For strike-slip faults 5.94±0.06 ≤ M _w ≤ 6.76±0.01, while for normal slip faults 5.78±0.13 ≤ M _w ≤ 6.75±0.03. Uncertainties in calculated fault areas and fault widths are generally smaller for strike-slip faulting compared to normal slip faulting. Except for the shortest SRL considered, rupture width are greater for strike-slip kinematics compared to normal slip.

Table 6 reports the stress drops along the AKFS for given earthquake moment magnitudes (Eq. 12). As expected, the magnitude of the stress drop does not change significantly with increasing moment magnitude and we derive an average stress drop of 2.6±0.2 MPa for the range of magnitudes considered.

Coulomb failure stress change accumulation due to gradual slip of the AKFS as function of time since its most recent earthquake (5.450 ka±0.16 BP) is summarised in Figure 3. We calculate a stress accumulation rate of 450±10 Pa yr⁻¹ for the lower-bound estimate of the right-lateral slip rate (0.58 mm yr⁻¹). In this case a ΔCFS (equivalent to the stress drop calculated in Eq. 13 and reported in Table 6) of 2.6±0.2 MPa is reached nowadays. The calculated stress accumulation for the upper-bound value of the right-lateral slip (3.35 mm yr⁻¹) is 2700±67 Pa yr⁻¹. As a result, the ΔCFS for fault rupture would have been attained ∼1 ka after the most recent earthquake and thus around ∼4.5 ka ago. The ΔCFS for pure normal rupture (Table 6) would have been reached at the greater rate of strike-slip within 500 years after the most recent earthquake, and thus ∼2.5 ka ago at the smaller rate of strike-slip. The accuracy and precision of these results are subject to the accuracy of geochronological and palaeoseismological data reported by Kürçer and Gökten (2012) to estimate the timing of the most recent earthquake on the fault segment. Between 2.5±0.2 and 15±0.5 MPa of stress has accumulated on the fault since then, if both the timing and the range of right-lateral slip rates are taken as true.

Coulomb failure stress changes on fractures and hence potential magma pathways above the magma reservoir are summarised in Figures 4A–D. gradual slip on the fault segment in general causes unclamping of fractures perpendicular to the strike of the fault plane, but clamping parallel to it. The magnitude of ΔCFS has been increasing as a function of time since the last earthquake for both slip distributions and the magnitude of clamping of fractures is fractionally greater for the smaller dextral gradual slip rate (<2.5% difference for all magnitudes of future earthquakes). For fractures perpendicular to the strike of the fault, the rate of unclamping is 4.0 × 10⁻⁴ MPa yr⁻¹ for the smaller rate of right-lateral slip, and 2.3 × 10⁻³ MPa yr⁻¹ for the greater. Clamping of fractures parallel to the strike of the fault plane occurs at a rate of 3.1 × 10⁻⁴ MPa yr⁻¹ for the smaller rate of right-lateral slip, and 2.1 × 10⁻³ MPa yr⁻¹ for the greater.

Coulomb failure stress changes associated with the hypothetical fault ruptures exacerbate the trend of unclamping of fractures perpendicular to the strike of the fault, and clamping of fractures parallel to it. Again, there is little difference in results between the two rates of strike-slip kinematics (<0.05 MPa). The magnitude of ΔCFS largely increases with magnitude of the hypothetical earthquake event. An earthquake with an SRL of 10 km (M _w ∼6.28) anomalously produces a smaller degree of clamping on fractures parallel to the strike of the fault plane than an earthquake with SRL of 5 km (M _w ∼5.94).

Earthquakes with magnitudes 6.48 ≤ M _w ≤ 6.76 significantly exacerbate the trends predicted by gradual Coulomb failure stress changes before fault rupture, whereby fractures perpendicular to the AKFS continue to be unclamped while fractures parallel to the AKFS being further clamped.

Figures 4E–H depict the Coulomb failure stress changes on fractures below Hasan Dağ’s magma reservoir due to gradual slip along the fault segment. Again, fractures are in general unclamped in a direction perpendicular to the fault plane but clamped parallel to it. There is <0.05 MPa difference in accrued stresses for the two slip rates for corresponding SRLs. Unclamping of fractures below the reservoir follows the same temporal pattern as predicted for fractures above the reservoir. The smaller slip-rate unclamps fractures perpendicular to the fault plane at a rate of 3.9 × 10⁻⁴ MPa yr⁻¹, while the greater slip-rate unclamps at a rate of 2.2 × 10⁻³ MPa yr⁻¹. Clamping of fractures parallel to the strike of the fault plane and below the reservoir occurs at a rate of 9.2 × 10⁻⁵ MPa yr⁻¹ for the smaller component of right-lateral slip, and 1.6 × 10⁻³ MPa yr⁻¹ for the greater.

All hypothetical earthquakes produce positive values of Coulomb failure stress changes on fractures in a direction perpendicular to the fault and hence result in further unclamping of potential magma pathways. This magnitude of clamping increases with magnitude of earthquake. All hypothetical earthquakes also exacerbate the trend of clamping on the fractures parallel to the strike of the fault plane. Matching the trend for fractures above the chamber, earthquakes with an SRL of 10 km produce a smaller magnitude of Coulomb failure stress change than earthquakes with an SRL of 5 km. Earthquakes with a SRL greater than 10 km produce increasing changes in Coulomb failure stress with increasing magnitude. Variations in the amount of modeled dextral slip on the fault segment have limited effects on the spatial and temporal evolution of Coulomb failure stress change on magma pathways below the reservoir.

The two rates of gradual slip produce very similar magnitudes of clamping of fractures along the reservoir wall; however, they are between one and two orders of magnitude smaller than the Coulomb failure stress change of fractures above and below the chamber (Figure 5). Point probes c and d (i.e., with tangents parallel to the fault plane) are both predicted to undergo clamping during gradual slip and prior to all hypothesised earthquake events. Point probes e and f (i.e., with tangents perpendicular to the fault plane) undergo unclamping during gradual slip.

The pattern of Coulomb failure stress change produced on fractures parallel to the reservoir produced by the hypothetical earthquakes is more complex. While both slip rates produce very similar magnitudes of Coulomb failure stress change (<0.03 MPa), earthquakes with SRL = 5, 20, 25 and 27 km (M _w = 5.94, 6.62, 6.73 and 6.76) cause clamping at c, d, and f. An earthquakes with a SRL of 15 km (M _w = 6.48) clamps c and f, while an earthquake with SRL = 10 km (M _w = 6.62) only clamps f. The magnitude of Coulomb failure stress change as a result of hypothetical earthquakes is much greater than that due to gradual slip on the fault. Therefore the spatial pattern of stress changes produced by seismic slip of the fault dominates the cumulative Coulomb failure stress changes at each point probe on the reservoir.

Figure 6 shows the stress change on the KNSZ at different times (1, 4 and 7 ka since the most recent earthquake) and depths (4, 6 and 8 km) due to gradual slip on the AKFS. Due to its greater plausibility, we only used the smaller estimate of right-lateral slip on the AKFS in the calculations. Unclamping of the KNSZ is predicted to propagate from the intersection with the AKFS to the southwest as a function of time. To the southeast of the intersection, the shear zone is clamped. While the distributions of ΔCFS show little variation with depth, the magnitude of ΔCFS decreases with distance from the AKFS and Hasan Dağ’s magma reservoir.

Coulomb failure stress changes on the KNSZ and vertical receiver faults parallel to it as a result of hypothesised seismic slip on the AKFS are shown in Figure 7. Coulomb failure stress changes are shown for depths of 4, 6, and 8 km by three earthquakes with SRLs of 5, 15, and 25 km, respectively. Overall, magnitudes of Coulomb failure stress change increase as a function of SRL. At all depths considered, interactions with the pressurised magma reservoir distort the pattern of Coulomb failure stress change, though this interaction becomes less prevalent with increasing magnitude of the earthquake. As a result of an M5.94 earthquake and a SRL of 5 km, the shear zone is predominantly unclamped to the west of the intersection with the TGFZ at all depths. At 4 km depth stress interactions are most significant and while the shear zone is broadly clamped to the southwest of the intersection up to ∼20 km distance, localised unclamping is predicted in the vicinity of the magma reservoir. These complexities also persist at greater depth.

The spatial patterns of Coulomb change for earthquakes of magnitudes 6.48 and 6.73 with SRLs of 15 and 25 km, respectively, show large similarities to each other. Clamping of the shear zone occurs to the southeast of the intersection (for an SRL = 15 km also to the northeast), and unclamping occurs to the southwest.

Earthquake magnitudes, fault dimensions and surface rupture lengths of hypothetical earthquakes along the AKFS were calculated from the Wells and Coppersmith (1994) empirical formulae. That study used a substantial number of intraplate source parameters giving statistically significant relationships (more than 95% probability level) between most faulting characteristics. The faulting kinematics (normal, reverse, strike-slip) proves insignificant for the resultant moment magnitude of potential earthquakes along the AKFS ( < 5% difference in results). However, it does prove significant for rupture areas associated with surface rupture lengths (SRLs) of 20 km and less. Given that the right-lateral component of slip along the fault segment is dominant over the normal slip, the empirical formulae for strike-slip faulting were deemed suitable to parameterise the fault dimensions. Calculating fault dimensions with the assumption of pure normal slip produces consistently smaller values. Omission of the normal slip component is likely to have caused a slight overestimation of the fault dimensions attributed to hypothetical earthquakes on the AKFS. In all models, we assumed rupturing to occur to the surface, reflecting rupture conditions caused by previous earthquakes studied by Kürçer and Gökten (2012).

Changes in Coulomb failure stress are static changes as changes in the off-fault stress field predicted by our models due to rupture of the AKFS are modeled instantaneous and permanent. Static stress change amplitudes diminish with 1/R³, where R is the distance from the earthquake’s epicentre (Walter and Amelung, 2007). The proximity of the volcano’s magma reservoir to the AKFS implies that seismic slip will be a dominant driver of stress changes in the reservoir as also demonstrated by results shown in panels B and E of Figure 5.

The two values of gradual right-lateral slip rates on the AKFS explored in the models cover the range of published values Kürçer and Gökten (2012). For the upper-bound estimate of right-lateral displacement (3.35 mm yr⁻¹) the models predict a seismic event of any plausible magnitude should have already occurred on the AKFS. If all model parameters are taken as exact representations of the geodynamic setting around Hasan Dağ, then this right-lateral component of slip for the AKFS must be regarded an overestimation. Our analysis indicates that the smaller component of right-lateral slip (0.58 mm yr⁻¹) likely provides more plausible estimate of the fault’s kinematics. In this case stress accumulation on the fault plane since the most recent earthquake (∼5.5 ka B.P.) is at least 2.5± 0.20 MPa and sufficient to promote failure of the segment (Figure 3).

Kürçer and Gökten (2012) propose a recurrence interval of 4660 ± 410 years for the AKFS. A greater quantity of time has elapsed since the last earthquake, even accounting for the upper boundary of error. Such a disparity in results may be due to a sparse data set available to the authors, but also due the authors’ assumption that the fault segment is dominated by normal slip with a minor right-lateral component. This contrasts the assumptions of fault kinematics behind our study based on the findings of Aktuğ et al. (2013) and Krystopowicz et al. (2020). Furthermore the discrepancy in findings can be explained by the evolution of Hasan Dağ’s magmatic system since its most recent eruption ∼9 ka B.P (Schmitt et al., 2014). With the longevity of the Hasan Dağ magmatic system, it is not unreasonable to assume steady-state replenishment of the magma reservoir matching that of the Holocene eruption rate predicted by Friedrichs et al. (2020) (0.1 km³ ka⁻¹). Pressurisation of the magma reservoir over this period would increase the normal stress on the AKFS, and thus discouraging slip.

In our model, Hasan Dağ’s magma reservoir is kept at constant (lithostatic) pressure (i.e., no pressure changes over time are invoked. Exact quantification of the current state of pressurisation of Hasan Dağ’s magma reservoir is beyond the scope of this study and possible variations in reservoir pressure due to magma influx and rejuvenation or volatile loss are hence neglected. Although InSAR data from the European Space Agency Sentinel-1 mission between 2014 and 2020 show annual ground displacement velocities at Hasan Dağ below the detectabilty limit of 1 cm yr⁻¹ (Biggs et al., 2021), the volcano’s high prominence and crustal mechanics of the CAVP could contribute to significant volumes of stealth magma accumulation at depth (Males and Gottsmann, 2021) at the geodetic detectability limit.

Hasan Dağ’s magma reservoir is sufficiently close to the fault segment to permit stress interactions between the two. Globally, there is statistical significance between earthquakes of M ≥ 5 and near field (within a few tens of km) volcanic unrest (Lemarchand and Grasso, 2007). Concurrency between earthquakes and eruptions is very high if M > 6 and eruptions are voluminous (Walter and Amelung, 2006). All point probes on the magma reservoir wall show linear increases, or decreases, of Coulomb failure stress as a result of protracted gradual slip on the fault segment before the occurrence of any hypothetical earthquake. All modeled earthquakes cause a continuation of these trends at the point probes, either exacerbating conditions for reservoir failure and potential eruption or dampening them. Prior to all hypothetical earthquakes and independent of slip rate, magma pathways above the reservoir are unclamped in a direction perpendicular to the fault plane. When only strike-slip movement is considered along the fault (with a trend of 34°NW) then the magma reservoir sits in the tensional stress quadrant since the reservoir is located towards the south-eastern edge of the fault segment and within the western fault block where faulting would be encouraged (see Figures 7, 8 and also Section 5.4). Clamping on the same magma pathways occurs parallel to the fault plane due to the normal slip component of faulting.

Exacerbation of conditions allowing for magma propagation towards the surface does, however, not guarantee a volcanic eruption. Gudmundsson and Philipp (2006) concluded that for magma propagation to reach the surface the stress field around the pathways must be homogenised and favour extension fracture formation. At Hasan Dağ, the stability of tensional stresses parallel to the AKFS prior to all hypothetical earthquakes is favourable to magma propagation towards the surface, thus promoting conditions for eruption. The magnitude of gradual unclamping perpendicular to the AKFS is greater than gradual clamping orthogonal to it. The earthquake magnitudes considered here predicted to cause further unclamping of pathways above the reservoir and perpendicular to the AKFS thus continuing the encouragement of a future eruption. If the reservoir were close to critical failure, small changes in stress (±0.01 MPa) could trigger an eruption (Nostro et al., 1998; Kriswati et al., 2019). The magnitude of unclamping by the modeled earthquakes is predicated to be > 0.13 MPa and therefore significant.

Gradual slip along the fault segment also causes unclamping of fractures below the reservoir. All hypothetical earthquakes further exacerbate this, irrespective of the magnitude of right-lateral slip with unclamping magnitudes > 0.115 MPa and therefore significant. Unclamping of fractures beneath Hasan Dağ’s magma reservoir (predicted to occur perpendicular to the AKFS) may already be evidenced by earthquakes with hypocentral depths > 10 km and thus below the average penetration depth of the TGFZ) (Kürçer and Gökten, 2014). Dilation of fractures below the reservoir may facilitate the formation of feeder dykes from depth and may encourage reservoir failure and eruption by the thermal and volumetric effects of injecting deeper-sourced mafic magma into an upper-crustal reservoir (Sparks et al., 1977). Mineralogical variations in ignimbrites at Hasan Dağ evidence such magma mixing prior to previous explosive eruptions (Aydar and Gourgaud, 1998).

Calculated Coulomb failure stress changes on the reservoir walls by gradual slip of the AKFS show a systematic pattern of clamping of fractures perpendicular to the AKFS and unclamping parallel to it. However, the spatial pattern of Coulomb failure stress change due to future hypothetical earthquakes does not show discernible trends at the lateral point probes. No model predicts wholesale encouragement or discouragement of lateral reservoir failure from predicted stress changes accompanying fault rupture. Lateral expulsion of magma from the reservoir is plausible if unclamping of the reservoir walls overcome the tensile strength of the country rock. Sill injection into country rock reduces the tensile stresses within the vicinity of the reservoir (Gudmundsson, 2000) moving the system further from critical failure. Direct squeezing of the reservoir is debated as to its actual effectiveness at encouraging eruptions due to the possible inadequacy of the volume of magma mobilised during clamping (Eggert and Walter, 2009).

Coulomb failure stress changes on the reservoir walls predicted by our models are a function of static stress transfer. This is not a true reflection of the stress perturbations occurring in the near field immediately after an earthquake. Earthquakes are accompanied by the release of energy in the form of seismic waves. The propagation of these waves causes short-lived, yet high amplitude, dynamic stress changes (Walter and Amelung, 2007). In addition to the static Coulomb failure stress changes that might encourage eruption, dynamic stress transfer can exacerbate the conditions for an eruption by causing expansion of the gaseous phase within a reservoir (Walter et al., 2007) or wholesale disruption of a magma mush (Davis et al., 2007; Gottsmann et al., 2009). Walter et al. (2007) concluded that dynamic stress transfer was at least as critical as static Coulomb failure stress transfer in triggering the Mt. Merapi eruption in 2006. Consequences of dynamic stress transfer for Hasan Dağ’s magma reservoir are a topic for further research.

On the September 20, 2020, a M5.1 earthquake located 35 km to the southwest of Hasan Dağ’s summit at a depth of ∼8 km reactivated a segment of the previously thought inactive Karaman-Niğde shear zone (KNSZ) (Aydar et al., 2021). The earthquake and subsequent aftershocks occurred in a volume predicted by this study to have undergone unclamping due to gradual slip along the AKFS (Figure 8). This supports the results of our modelling.

Our study does not model changes in Coulomb failure stress on the magma reservoir and fault segment due to slip events along the KNSZ. Further exploration is needed to quantify the mid to long-term effects of these earthquakes. An earthquake occurring on the AKFS would cause rotations of the stress field encompassing the KNSZ and would therefore encourage, or discourage further earthquakes there. However, the orientation of the shear zone relative to the fault segment’s Coulomb failure stress change pattern means any seismic event along that segment only transfers a significant magnitude of Coulomb failure stress up to a distance ∼20–30 km from the intercept of the two features.

An empirical formula (Wang et al., 2006) links earthquake magnitude with maximum distance of liquefaction R_max log R max = 2.05 + 0.45 M w (14)

Taking the largest possible earthquake that could occur along the AKFS (M _w ∼ 6.76), the equation predicts an R_max of ∼120 km. This radius encompasses the towns of Aksaray, Nevşehir and Niğde amongst other settlements. These cities alone have a combined population of over 1,000,000. Furthermore, over a third of the city of Aksaray has been deemed to have a moderate to very high potential for liquefaction (Yalcin et al., 2008). Detailed estimations on liquefaction potential for the rest of the region is lacking.

Hasan Dağ is not the only volcanic centre that could be affected by seismic activity along the AKFS. Delle Donne et al. (2010) link earthquake magnitude to the greatest distance seismic energy released by that earthquake can trigger volcanic unrest: M w = − 6.4 + 2.17 l o g R (15) where R is the maximum distance of unrest. For the greatest earthquake magnitude predicted for the fault segment (M _w = 6.76), R is ∼ 1160 km. This encompasses the Quaternary volcanic centres of Acigöl-Nevsehir, Göllü Dag, Erciyes Dağ, Karadağ and the Karapinar volcanic field.

At Hasan Dağ, future eruptions are likely to follow the historic pattern of volcanic activity. This could express itself in two ways. Lava dome extrusions as a continuation of recent activity (Aydar and Gourgaud, 1998; Schmitt et al., 2014; Friedrichs et al., 2020) or a fourth catastrophic ignimbrite eruption capable of reaching Aksaray (Aydar and Gourgaud, 1998). Late Pleistocene and Quaternary dome formations indicate eruptive volumes of up to 0.37 km³ (Friedrichs et al., 2020) with pyroclastic flows extending up to 9 km from Hasan Dağ’s summit (Aydar and Gourgaud, 1998).

In summary, significant seismic and volcano-tectonic hazards exist along the AKFS and we recommend local risk assessment include scenarios explored in this study.