Investigating Stress Transfer Between the Tuz Gölü Fault Zone and Hasan Dağ Volcano (Turkey)

Faulting, magmatism and volcanism are intrinsically linked by plate tectonics. Fault slip imparts stress changes to the surrounding crust and other faults and fractures. Volcano-tectonic hazard assessment in areas with long recurrence intervals of volcanic and tectonic activity requires an assessment of current stress levels. Here we investigate stress transfer between the Akhisar-Kiliç fault segment (AKFS) of the seismically active Tuz Gölü fault zone in the Central Anatolian Volcanic Province and the active Hasan Dağ volcanic complex. Current stress accumulation by protracted gradual slip on the AKFS since its most recent rupture (5.45 ka±0.16 BP) is quantified using the Coulomb Failure Stress change (ΔCFS). We calculate currently accrued ΔCFS magnitudes of between 2.5±0.2 and 15±0.5 MPa on the fault plane for published lower and upper-bound estimates of right-lateral slip rates, respectively. These changes are sufficient to promote failure of the segment. The M5.1 September 20, 2020 earthquake SW of Hasan Dağ occurred in a volume predicted by this study to have undergone fault unclamping by gradual slip along the AKFS. We also show that gradual slip of the AKFS contributes to the progressive unclamping of fractures and transtensional opening of potential magma pathways oriented perpendicularly to the AKFS both above and below Hasan Dağ’s magma reservoir, while pathways oriented parallel to the AKFS are being clamped. Earthquake moment magnitudes of between M_w5.94 and M_w6.76 due to hypothetical partial or complete rupture of the segment exacerbate these trends. The spatial pattern of Coulomb failure stress changes on Hasan Dağ’s magma reservoir is predominantly controlled by the location of rupture relative to the reservoir with the magnitude of the earthquake playing a subordinate role. We explore implications of our findings for the assessment of interconnected seismo-volcanic hazards and associated risks.

1 Introduction

At local scales, the temporal and spatial patterns of volcanic and seismic events are stochastic, statistically unrelated to one another (Hill et al., 2002). However, faulting and magmatic processes are intrinsically linked by plate tectonics (e.g., Perfit and Davidson, 2000; Karaoğlu et al., 2017) and there are several examples of stress transfer between tectonic and magmatic systems where fault rupture led to volcanic activity (Table 1). The inflation or voiding of magma reservoirs before and after eruptions causes local stresses to rotate, bringing optimally orientated faults in areas of increased stress closer to failure (Fukuyama et al., 2001; Roman et al., 2004; Ruppert et al., 2011). Conversely, stress released via slip on faults can encourage failure of a magma reservoir (Gudmundsson and Brenner, 2003; Walter and Amelung, 2007; Walter et al., 2007; Bonali et al., 2013). A positive feedback loop between the two phenomena has been noted by Nostro et al. (1998). Often both phenomena are only linked in hindsight.

TABLE 1

TABLE 1. Selected examples of Coulomb failure stress transfer to magmatic systems causing the triggering (or partial triggering) of eruptions. Three modes are proposed: 1) unclamping of magma conduits above reservoir, 2) unclamping of magma conduits below reservoir, 3) direct clamping of the magma reservoir due to displacement on fault planes. Moment magnitude of earthquake (M_w), horizontal distance (D) distance between earthquake epicenter and volcanic system.

Reservoir failure and transport of magma from a reservoir to the surface are prerequisites for a volcanic eruption. Reservoir failure occurs when excess magmatic pressure above lithostatic exceeds the surrounding rock’s tensile strength (Wilson, 1980; Gudmundsson, 1988). Changes in magma reservoir volume typically lead to pressure variations on the reservoir wall (Blake, 1981). Equally, stress changes on a magma reservoir may produce volumetric fluctuations if the magnitude of stress change is great enough (Stein, 1999). Stress variations can lead to the clamping or unclamping of fractures around a reservoir which may promote or inhibit the propagation of magma towards the surface. Here, we investigate stress transfer between a major strike-slip fault system (the Tuz Gölü fault zone) and the active Hasan Dağ volcanic complex of the Central Anatolian Volcanic Province. A seismic swarm culminating in a M5.1 to M5.2 earthquake on September 20, 2020 occurred within 10 km of the summit of Hasan Dağ and within 20 km from the Tuz Gölü fault zone, in a volume with no recorded historic seismicity (Aydar et al., 2021). This raises questions about possible stress transfer from the major fault zone to the surrounding crust and secondary fault systems. We show that the currently accrued stress along the Tuz Gölü fault zone forces the magma-tectonic system towards failure and discuss resulting hazards.

2 Background

2.1 Central Anatolian Tectonics

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).

FIGURE 1

FIGURE 1. Maps of the study area. Inset: Map of tectonic setting of the Anatolian microplate bound by the Eurasian plate to the north and African and Arabian plates to the south. Large arrows next to plate labels show the overall direction of plate movement. Black lines represent fault zones (EFZ: Ecemiş fault zone; TGFZ: Tuz Gölü fault zone). The red triangle shows the location of Hasan Dağ volcano. Fault locations and plate directions are from Tosun et al. (2007). Main panel: solid red line: Akhisar-Kiliç fault segment of the TGFZ; dashed red line: approximation of the Karaman-Niğde shear zone; BHD: Big Hasan Dağ; SHD: Small Hasan Dağ.

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).

2.2 Hasan Dağ Volcano

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).

2.3 Volcano-Tectonics of Hasan Dağ

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).

3 Methods

Slip on a fault occurs when shear stress overcomes static friction across the fault plane (Wang et al., 2014). This principle underlies the Coulomb Failure Criterion, which predicts the occurrence of slip when the Coulomb failure stress change (ΔCFS) on a fault overcomes a given value, defined as (King et al., 1994):

$\mathrm{\Delta }CFS=\mathrm{\Delta }\tau +\mu \left(\mathrm{\Delta }{\sigma }_n+\mathrm{\Delta }P\right)(1)$

where Δτ is the change in shear stress, μ is the coefficient of friction, Δσ_n is the change in normal stress, and Δ P is the change in pore fluid pressure. Symbols and abbreviations used throughout this paper are summarised in Table 2. The conventions employed in this study are 1) shear stress change on a fault is positive in the direction of slip and 2) normal stress change is extensional and thus positive if a fault or fracture is unclamped. Failure is encouraged if ΔCFS is positive and discouraged if negative; both increased shear and unclamping of faults promote failure.

TABLE 2

TABLE 2. Abbreviations, units and definitions.

Assuming rocks are undrained during the coseismic phase the change in pore pressure is defined by (Cocco and Rice, 2002):

$\mathrm{\Delta }{\sigma }_p=-B\left(\mathrm{\Delta }{\sigma }_{kk}\right)/3(2)$

where B is Skempton’s coefficient and σ_kk is the trace of the stress tensor (Harris, 1998). Assuming that the ductility of the interior fault material is greater than that of the surrounding materials gives σ_xx = σ_yy = σ_zz, and (Δσ_kk)/3 = Δσ (Harris, 1998). These assumptions cause the simplification of Eq. 1 to:

$\mathrm{\Delta }CFS=\mathrm{\Delta }\tau -{\mu }^{\prime }\mathrm{\Delta }\sigma (3)$

where is the effective coefficient of friction (King et al., 1994):

${\mu }^{\prime }=\mu \left(1-B\right)(4)$

In this study μ′ = 0.4, in line with volcano-tectonic studies of similar sized fault systems (Nostro et al., 1998).

3.1 Model Set-Up

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/).

FIGURE 2

FIGURE 2. Set up of the 3D finite element model. (A) Layered crustal mechanical heterogeneities are represented by subdomains with properties defined in Table 3. The modelled subsurface architecture shows an increase in mechanical rigidity and density with depth. TGFZ: Tuz Gölü fault zone; KNSZ: Karaman-Niğde shear zone. An oblate cavity [shown in (B)] centred at [0, 0, −5000], but not shown for clarity represents the magma reservoir of Hasan Dağ. (B) Schematic (not to scale) of the modeled relationship between the AKFS’s fault plane and Hasan Dağ’s magma reservoir. The magma reservoir is located towards the south-eastern end of the fault plane.

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+blog\left(SRL\right)(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:

$log\left(RA\right)=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:

$log\left(RW\right)=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.

3.2 Benchmarking

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.

3.3 Crustal Heterogeneity

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):

$\rho =1.6612V_p-0.4721V_p^2+0.0671V_p^3-0.0043V_p^4+0.000106V_p^5(8)$

Poisson’s ratio (ν) was calculated as a function of the p- and s-wave velocity ratios:

$\nu =\left({\left[V_p/V_s\right]}^2-2\right)/\left({\left[V_p/V_s\right]}^2-1\right)(9)$

Finally, the dynamic Young’s modulus (E) was calculated from:

$E_d=\left(V_p^2\nu \left(1+\nu \right)\left(1-2\nu \right)\right)/\left(1-\nu \right)(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.

TABLE 3

TABLE 3. Mechanical properties assigned to each subdomain in the model. The plateau surrounding Hasan Dağ is set to z = 0 m in the model, thus the edifice extends to a height of z = 2532 m. See Eq. 8 and Supplementary Table S1 for further information on the derivation of the values.

3.4 Stress Drop of Earthquakes

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\mathit{log}\left(M_o\right)-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):

$\mathrm{\Delta }\sigma =\left(2M_o\right)/\left(\pi R{W}^{2}SRL\right)(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.

3.5 Time-Dependent Model

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.

3.6 Coulomb Failure Stress Change on Receiver Faults

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.

3.7 Average Displacement by Faulting

The average fault displacement AD due to each earthquake was calculated via

$log\left(AD\right)=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.

TABLE 4

TABLE 4. Values of average displacement (AD) (1σ errors in parentheses) across the Akhisar-Kiliç fault segment for given surface rupture length (SRL) derived from Eq. 13. Strike-slip and normal slip components are calculated based on the ratio between dextral slip and normal slip components of gradual annual slip on the segment over the past 2 Ma for given lower 1) and upper 2) bound slip rates as reported in Krystopowicz et al. (2020).

3.8 Coulomb Failure Stress Change on Hasan Dağ’s Magma Reservoir

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.

4 Results

4.1 Fault Parameters

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 5

TABLE 5. Values of maximum moment magnitude M_w for each earthquake for selected surface rupture lengths (SRL), the associated rupture areas (RA) and rupture widths (RW) as given by Eqs 5, 6. 1σ errors are reported in parentheses. * gives values for fault segment with pure strike-slip kinematics; ** gives values for fault segment with pure normal slip kinematics.

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.

TABLE 6

TABLE 6. Values of stress drop Δσ corresponding the moment magnitudes M_w of earthquakes from Table 5. * for fault segment with pure strike-slip kinematics. ** for fault segment with pure normal slip kinematics. The seismic moment M_o was derived from the moment magnitude M_w of each earthquake via Eq. 11 (Hanks and Kanamori, 1979) which gives seismic moment in dyn cm, converted to MPa (1 dyn cm = 1 × 10⁻⁷ MPa) and substituted into Eq. 12. 1σ errors are reported in parentheses.

4.2 Coulomb Failure Stress Changes on Akhisar-Kiliç Fault Segment

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.

FIGURE 3

FIGURE 3. Modeled Coulomb failure stress change accumulation ΔCFS (Eq. 3) by gradual slip on the Akhisar-Kiliçc fault segment (AKFS) for the range of transtensional slip velocities reported by Krystopowicz et al. (2020) as a function of time elapsed since the most recent earthquake [5450 ± 160; using data from Kürçer and Gökten (2012)]. The vertical shaded area represents modern-day including data uncertainty and the horizontal shaded areas indicates stress drops Δσ of 2.6±0.2 MPa and 1.2±0.2 MPa calculated for the range of moment magnitudes M_w of earthquakes reported in Table 5. The solid blue line shows the stress accumulation on the AKFS for a dextral slip rate of 0.58 mm yr⁻¹ (labeled slow rate) with a gradient of 450±10 Pa yr⁻¹. The solid black line shows the stress accumulation for a dextral slip rate of 3.35 mm yr⁻¹ (labeled fast rate) and a resultant gradient of 2700 ± 67 Pa yr⁻¹. Both scenarios have a normal slip rate of 0.05 mm yr⁻¹. Uncertainties are expressed by dotted lines.

4.3 Coulomb Failure Stress Change on Hasan Dağ’s Magma Reservoir and Feeder Systems

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.

FIGURE 4

FIGURE 4. Bar graphs of modeled changes in Coulomb failure stress on vertical fractures above (A–D) and below (E–H) Hasan Dağ’s magma reservoir (extending from and to point probes a and b, respectively, in Figure 2B). The blue bars represent the change in Coulomb failure stress due to gradual slip along the Akhisar-Kiliçc fault segment (AKFS) since its most recent rupture (5450 ± 160 ka) and prior to hypothetical earthquakes with given moment magnitude (M_w). Orange segments show the Coulomb failure stress change as a result of the hypothetical earthquakes. (A,B,E,F) show results for a rate of 0.58 mm yr⁻¹ right-lateral slip and a rate of 0.053 6 mm yr⁻¹ of normal slip on the fault segment (labeled slow slip), while (C,D,G,H) show Coulomb failure stress change due to right-lateral slip of 3.35 mm y⁻¹ and 0.053 6 mm yr⁻¹ of normal slip (labeled fast slip). Positive stress changes represent unclamping of fractures, while negative stress changes indicate clamping. (A,C,E,G) show stress changes perpendicular to the strike of the AKFS. (B,D,F,H) show stress changes in a direction that is parallel to the strike of the AKFS.

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.

FIGURE 5

FIGURE 5. Same as Figure 4 but for modeled Coulomb failure stress changes on point probes c, d, e and f on Hasan Dağ’s magma reservoir (see also Figure 2B). Tangents of c and d are parallel to the AKFS while tangents on e and f are perpendicular to it. (A,D) show the Coulomb failure stress changes on the points due to the gradual stress accumulation prior to hypothetical earthquakes with given moment magnitude. (B,E) show the Coulomb failure stress changes as a result of hypothetical earthquakes with given SRL. (C,F) show the cumulative Coulomb failure stress changes. (A–C) for 0.58 mm yr⁻¹ right-lateral slip and 0.053 6 mm yr⁻¹ of normal slip (labeled slow slip). (D–F) for 3.35 mm yr⁻¹ right-lateral slip and 0.053 6 mm yr⁻¹ of normal slip (labeled fast 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.

4.4 Coulomb Failure Stress Changes on Receiver Faults of the KNSZ

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.

FIGURE 6

FIGURE 6. Modeled gradual Coulomb failure stress changes at selected depths on vertical planes parallel to the Karaman-Niğde shear zone (trending 55° NE) after 1000, 4000 and 7000 years since the most recent earthquake on the AKFS. Changes are calculated for the more plausible slow gradual slip rate of 0.58 mm yr⁻¹ right-lateral and 0.053 6 mm yr⁻¹ normal slip. Coulomb failure stress changes at 4 km depth (A–C), 6 km depth (D–F) and 8 km depth (G–I). Note the reduction in scale in panel (G).

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.

FIGURE 7

FIGURE 7. Modeled Coulomb failure stress changes at selected depths on vertical planes parallel to the Karaman-Niğde shear zone (trending 55° NE) as a result of hypothetical earthquakes along the AKFS with a rake of −175° commensurate with the gradual slip distribution of 0.58 mm yr⁻¹ right-lateral slip, and 0.053 6 mm yr⁻¹ normal slip. (A–C) Coulomb failure stress change due to aM_w = 5.94 earthquake with surface rupture length (SRL) of 5 km. (D–F) for M_w = 6.48 and SRL = 15 km. (G–I) for M_w = 6.73 and SRL = 25 km. Note the change in scales as magnitude and depth increase.

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.

5 Discussion

5.1 Fault Parameters

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.

5.2 Rupture of the Akhisar-Kiliç Fault Segment

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.

5.3 Coulomb Failure Stress Change on Hasan Dağ’s Magma Reservoir

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.

FIGURE 8

FIGURE 8. Modeled present-day accumulation of Coulomb failure stress change on vertical faults parallel to the Karaman-Niğde shear zone (KNSZ) at a depth of 8 km due to right-lateral gradual slip of 0.58 mm yr⁻¹ along the TGFZ. The black star marks the epicentral location (from the USGS earthquake catalogue) of the 20th of September 2020 M5.1 earthquake. The earthquake occurred in an volume of predicted positive Coulomb failure stress change where the shear zone would have undergone unclamping, thus encouraging rupture.

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.

5.4 The September 20, 2020 Karaman-Niğde Earthquake

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.

5.5 Implications for Natural Hazard and Risk Assessment

An empirical formula (Wang et al., 2006) links earthquake magnitude with maximum distance of liquefaction R_max

$\mathit{log}\left(R_{\max }\right)=2.05+0.45M_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.17log\left(R\right)(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.

6 Conclusion

Using 3D finite element modelling, this study shows that protracted gradual slip along the Akhisar-Kiliç fault segment (AKFS) of the Tuz Gölü fault zone forces the fault segment and Hasan Dağ’s magma reservoir towards failure. The current accrual of Coulomb failure stress changes along the fault segment progressively unclamps magma pathways both above and below the magma reservoir. For the upper-bound estimate of Coulomb failure stress change accumulation of ∼ 2700±67 Pa yr⁻¹ strike-slip failure conditions were reached ∼4.5 ka B.P. For the lower-bound estimate of 450±10 Pa yr⁻¹ the conditions are currently met for strike-slip failure and were met ∼2.5 ka ago for normal slip failure. Our findings imply that the lower-bound of published right-lateral slip rates (0.58 mm yr⁻¹) (Krystopowicz et al., 2020) might be a more plausible estimate of the slip distribution and stress accumulation along the fault segment. Its rupture could result in an earthquake of 5.78±0.13 ≤ M_w ≤ 6.76±0.03 in the future.

Notwithstanding fault rupture, the current stress transfer from the AKFS to Hasan Dağ’s magma reservoir by gradual slip of the segment alone is moving the reservoir towards failure with potential for renewed eruptive activity.

The criticality of Hasan Dağ’s magmatic system is unknown, but fumarolic degassing and mid-crustal seismicity anomalies are suggestive of active magmatic feeder systems beneath the volcano. These observations in combination with our findings highlight a multi-hazard scenario that is ideally considered in interconnected natural risk analysis and disaster risk reduction (Pescaroli and Alexander, 2018). Finally, our findings have implications for quantifying volcano-tectonic interactions in systems with a long paucity in tectonic and eruptive activity.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://doi.org/10.5281/zenodo.5036587.

Author Contributions

JG conceived and coordinated the research underpinning this paper. EH developed the models, analyses and initial interpretation as part of an undergraduate project and thesis at the University of Bristol under the supervision of JG. Both authors drafted, read and approved the final manuscript.

Funding

JG acknowledges support from the Newton Fund TurkVolc project, NE/P008437/1 and NERC grant NE/S008845/1.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

JG expresses his gratitude to TurkVolc project team members Gökhan Atıcı and Mehmet Çobankaya for stimulating discussions whilst climbing Hasan Dağ volcano that inspired the study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2021.732696/full#supplementary-material

References

Abgarmi, B., Delph, J. R., Ozacar, A. A., Beck, S. L., Zandt, G., Sandvol, E., et al. (2017). Structure of the Crust and African Slab beneath the Central Anatolian Plateau from Receiver Functions: New Insights on Isostatic Compensation and Slab Dynamics. Geosphere 13, 1774–1787. doi:10.1130/GES01509.1

CrossRef Full Text | Google Scholar

Aktuğ, B., Parmaksız, E., Kurt, M., Lenk, O., Kılıçoğlu, A., Gürdal, M. A., et al. (2013). Deformation of Central Anatolia: Gps Implications. J. Geodyn. 67, 78–96. doi:10.1016/j.jog.2012.05.008

CrossRef Full Text | Google Scholar

Aydar, E., and Gourgaud, A. (1998). The Geology of Mount Hasan Stratovolcano, Central Anatolia, Turkey. J. Volcanol. Geotherm. Res. 85, 129–152. doi:10.1016/S0377-0273(98)00053-5

CrossRef Full Text | Google Scholar

Aydar, E., Diker, C., Ulusoy, İ., and Şen, E. (2021). Volcanic Unrest Possibilities in Response to Recent Obruk Seismic Swarm on and Around Hasandağ Stratovolcano (Central Anatolia, Turkey). Comptes RendusGéosci. 353, 1–18. doi:10.5802/crgeos.46

CrossRef Full Text | Google Scholar

Aydar, E. (1992). Etude volcano-structurale et magmatologique du strato-volcan Hasan Dagi (Anatolie Centrale-Turquie). Clermont-Ferrand, France: Université Blaise Pascal. Thesis.

Google Scholar

Biggs, J., Dogru, F., Dagliyar, A., Albino, F., Yip, S., Brown, S., et al. (2021). Baseline Monitoring of Volcanic Regions with Little Recent Activity: Application of Sentinel-1 InSAR to Turkish Volcanoes. J. Appl. Volcanol. 10, 1–14. doi:10.1186/s13617-021-00102-x

CrossRef Full Text | Google Scholar

Blake, S. (1981). Volcanism and the Dynamics of Open Magma chambers. Nature 289, 783–785. doi:10.1038/289783a0

CrossRef Full Text | Google Scholar

Bonali, F. L., Tibaldi, A., Corazzato, C., Tormey, D. R., and Lara, L. E. (2013). Quantifying the Effect of Large Earthquakes in Promoting Eruptions Due to Stress Changes on Magma Pathway: the Chile Case. Tectonophysics 583, 54–67. doi:10.1016/j.tecto.2012.10.025

CrossRef Full Text | Google Scholar

Bonali, F. L. (2013). Earthquake-induced Static Stress Change on Magma Pathway in Promoting the 2012 Copahue Eruption. Tectonophysics 608, 127–137. doi:10.1016/j.tecto.2013.10.006

CrossRef Full Text | Google Scholar

Brocher, T. M. (2005). Empirical Relations between Elastic Wavespeeds and Density in the Earth's Crust. Bull. Seismol. Soc. Am. 95, 2081–2092. doi:10.1785/0120050077

CrossRef Full Text | Google Scholar

Cocco, M., and Rice, J. R. (2002). Pore Pressure and Poroelasticity Effects in Coulomb Stress Analysis of Earthquake Interactions. J. Geophys. Res. 107, ESE 2–1–ESE 2–17. doi:10.1029/2000JB000138

CrossRef Full Text | Google Scholar

Davis, M., Koenders, M. A., and Petford, N. (2007). Vibro-agitation of Chambered Magma. J. Volcanol. Geotherm. Res. 167, 24–36. doi:10.1016/j.jvolgeores.2007.07.012

CrossRef Full Text | Google Scholar

De la Cruz-Reyna, S., Tárraga, M., Ortiz, R., and Martínez-Bringas, A. (2010). Tectonic Earthquakes Triggering Volcanic Seismicity and Eruptions. Case Studies at Tungurahua and Popocatépetl Volcanoes. J. Volcanol. Geotherm. Res. 193, 37–48. doi:10.1016/j.jvolgeores.2010.03.005

CrossRef Full Text | Google Scholar

Delle Donne, D., Harris, A. J. L., Ripepe, M., and Wright, R. (2010). Earthquake-induced thermal Anomalies at Active Volcanoes. Geology 38, 771–774. doi:10.1130/G30984.1

CrossRef Full Text | Google Scholar

Deniel, C., Aydar, E., and Gourgaud, A. (1998). The Hasan Dagi Stratovolcano (Central Anatolia, Turkey): Evolution from Calc-Alkaline to Alkaline Magmatism in a Collision Zone. J. Volcanol. Geotherm. Res. 87, 275–302. doi:10.1016/S0377-0273(98)00097-3

CrossRef Full Text | Google Scholar

Dİker, C., Ulusoy, I., Akkaş, E., Şen, E., Çubukçu, H., Gümüş, E., et al. (2018). “Self-potential, Surface Temperature, CO₂ Measurements and Related Structural Discontinuities on Mount Hasan, Turkey,” in Cities On Volcanoes 10 (INGV).

Google Scholar

Eggert, S., and Walter, T. R. (2009). Volcanic Activity before and after Large Tectonic Earthquakes: Observations and Statistical Significance. Tectonophysics 471, 14–26. doi:10.1016/j.tecto.2008.10.003

CrossRef Full Text | Google Scholar

Fichtner, A., Saygin, E., Taymaz, T., Cupillard, P., Capdeville, Y., and Trampert, J. (2013). The Deep Structure of the North Anatolian Fault Zone. Earth Planet. Sci. Lett. 373, 109–117. doi:10.1016/j.epsl.2013.04.027

CrossRef Full Text | Google Scholar

Friedrichs, B., Atıcı, G., Danišík, M., Atakay, E., Çobankaya, M., Harvey, J. C., et al. (2020). Late Pleistocene Eruptive Recurrence in the post-collisional Mt. Hasan Stratovolcanic Complex (Central Anatolia) Revealed by Zircon Double-Dating. J. Volcanol. Geotherm. Res. 404, 107007. doi:10.1016/j.jvolgeores.2020.107007

CrossRef Full Text | Google Scholar

Fukuyama, E., Kubo, A., Kawai, H., and Nonomura, K. i. (2001). Seismic Remote Monitoring of Stress Field. Earth Planet. Sp 53, 1021–1026. doi:10.1186/BF03351699

CrossRef Full Text | Google Scholar

Göğüş, O. H., Pysklywec, R. N., Şengör, A. M. C., and Gün, E. (2017). Drip Tectonics and the Enigmatic Uplift of the Central Anatolian Plateau. Nat. Commun. 8, 1–9. doi:10.1038/s41467-017-01611-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Gottsmann, J., Lavallée, Y., Martí, J., and Aguirre-Díaz, G. (2009). Magma-tectonic Interaction and the Eruption of Silicic Batholiths. Earth Planet. Sci. Lett. 284, 426–434. doi:10.1016/J.Epsl.2009.05.008

CrossRef Full Text | Google Scholar

Gudmundsson, A., and Brenner, S. L. (2003). Loading of a Seismic Zone to Failure Deforms Nearby Volcanoes: a New Earthquake Precursor. Terra Nova 15, 187–193. doi:10.1046/j.1365-3121.2003.00481.x

CrossRef Full Text | Google Scholar

Gudmundsson, A., and Philipp, S. L. (2006). How Local Stress Fields Prevent Volcanic Eruptions. J. Volcanol. Geotherm. Res. 158, 257–268. doi:10.1016/j.jvolgeores.2006.06.005

CrossRef Full Text | Google Scholar

Gudmundsson, A. (1988). Effect of Tensile Stress Concentration Around Magma chambers on Intrusion and Extrusion Frequencies. J. Volcanol. Geotherm. Res. 35, 179–194. doi:10.1016/0377-0273(88)90015-7

CrossRef Full Text | Google Scholar

Gudmundsson, A. (2000). Dynamics of Volcanic Systems in Iceland: Example of Tectonism and Volcanism at Juxtaposed Hot Spot and Mid-ocean ridge Systems. Annu. Rev. Earth Planet. Sci. 28, 107–140. doi:10.1146/annurev.earth.28.1.107

CrossRef Full Text | Google Scholar

Hamling, I. J., and Kilgour, G. (2020). Goldilocks Conditions Required for Earthquakes to Trigger Basaltic Eruptions: Evidence from the 2015 Ambrym Eruption. Sci. Adv. 6, eaaz5261. doi:10.1126/sciadv.aaz5261

PubMed Abstract | CrossRef Full Text | Google Scholar

Hanks, T. C., and Kanamori, H. (1979). A Moment Magnitude Scale. J. Geophys. Res. 84, 2348–2350. doi:10.1029/JB084iB05p02348

CrossRef Full Text | Google Scholar

Harris, R. A. (1998). Introduction to Special Section: Stress Triggers, Stress Shadows, and Implications for Seismic hazard. J. Geophys. Res. 103, 24347–24358. doi:10.1029/98JB01576

CrossRef Full Text | Google Scholar

Hill, D. P., Pollitz, F., and Newhall, C. (2002). Earthquake-Volcano Interactions. Phys. Today 55, 41–47. doi:10.1063/1.1535006

CrossRef Full Text | Google Scholar

Kanamori, H., and Anderson, D. L. (1975). Theoretical Basis of Some Empirical Relations in Seismology. Bull. Seismol. Soc. Am. 65, 1073–1095.

Google Scholar

Karaoğlu, Ö., Selçuk, A. S., and Gudmundsson, A. (2017). Tectonic Controls on the Karlıova Triple junction (Turkey): Implications for Tectonic Inversion and the Initiation of Volcanism. Tectonophysics 694, 368–384. doi:10.1016/j.tecto.2016.11.018

CrossRef Full Text | Google Scholar

King, G. C., Stein, R. S., and Lin, J. (1994). Static Stress Changes and the Triggering of Earthquakes. Bull. Seismol. Soc. Am. 84, 935–953.

Google Scholar

Kiser, E., Palomeras, I., Levander, A., Zelt, C., Harder, S., Schmandt, B., et al. (2016). Magma Reservoirs from the Upper Crust to the Moho Inferred from High-Resolution Vp and vs Models beneath Mount St. Helens, Washington State, USA. Geology 44, 411–414. doi:10.1130/G37591.1

CrossRef Full Text | Google Scholar

Kürçer, A., and Gökten, Y. E. (2012). “Paleoseismological Three Dimensional Virtual Photography Method; A Case Study: Bağlarkayasi-2010 Trench, Tuz Gölü Fault Zone, Central Anatolia, Turkey,” in Tectonics-Recent Advances. Editor Sharkov, E. (Rijeka, Croatia: InTech,), 201–228.

Google Scholar

Kriswati, E., Meilano, I., Iguchi, M., Abidin, H. Z., and Surono, (2019). An Evaluation of the Possibility of Tectonic Triggering of the Sinabung Eruption. J. Volcanol. Geotherm. Res. 382, 224–232. doi:10.1016/j.jvolgeores.2018.04.031

CrossRef Full Text | Google Scholar

Krystopowicz, N. J., Schoenbohm, L. M., Rimando, J., Brocard, G., and Rojay, B. (2020). Tectonic Geomorphology and Plio-Quaternary Structural Evolution of the Tuzgölü Fault Zone, Turkey: Implications for Deformation in the interior of the Central Anatolian Plateau. Geosphere 16, 1107–1124. doi:10.1130/GES02175.1

CrossRef Full Text | Google Scholar

Kürçer, A., and Gökten, Y. E. (2014). Neotectonic-period Characteristics, Seismicity, Geometry and Segmentation of the Tuz Gölü Fault Zone. Bull.Min.Res.Exp. 149, 19–68. doi:10.19111/BMRE.94099

CrossRef Full Text | Google Scholar

Gencalioglu Kuscu, G., and Geneli, F. (2010). Review of post-collisional Volcanism in the Central Anatolian Volcanic Province (Turkey), with Special Reference to the Tepekoy Volcanic Complex. Int. J. Earth Sci. (Geol Rundsch) 99, 593–621. doi:10.1007/s00531-009-0504-710.1007/s00531-008-0402-4

CrossRef Full Text | Google Scholar

Lemarchand, N., and Grasso, J.-R. (2007). Interactions between Earthquakes and Volcano Activity. Geophys. Res. Lett. 34, L24303. doi:10.1029/2007gl031438

CrossRef Full Text | Google Scholar

Males, K., and Gottsmann, J. (2021). Minimum Detectable Mass and Volume Fluxes during Magmatic Recharge at High Prominence Volcanoes: An Application to Erciyes Dağ Volcano (Turkey). Front. Earth Sci. 9, 750063. doi:10.3389/feart.2021.750063

CrossRef Full Text | Google Scholar

Nostro, C., Stein, R. S., Cocco, M., Belardinelli, M. E., and Marzocchi, W. (1998). Two-way Coupling between Vesuvius Eruptions and Southern Apennine Earthquakes, italy, by Elastic Stress Transfer. J. Geophys. Res. 103, 24487–24504. doi:10.1029/98JB00902

CrossRef Full Text | Google Scholar

Perfit, M. R., and Davidson, J. P. (2000). “Plate Tectonics and Volcanism,” in Encyclopedia of Volcanoes. Editor Sigurdsson, H. San Diego, CA: Academic Press, 89–113.

Google Scholar

Pescaroli, G., and Alexander, D. (2018). Understanding Compound, Interconnected, Interacting, and Cascading Risks: A Holistic Framework. Risk Anal. 38, 2245–2257. doi:10.1111/risa.13128

PubMed Abstract | CrossRef Full Text | Google Scholar

Piromallo, C., and Morelli, A. (2003). P Wave Tomography of the Mantle under the Alpine-Mediterranean Area. J. Geophys. Res. 108. doi:10.1029/2002JB001757

CrossRef Full Text | Google Scholar

Roman, D. C., Power, J. A., Moran, S. C., Cashman, K. V., Doukas, M. P., Neal, C. A., et al. (2004). Evidence for dike Emplacement beneath Iliamna Volcano, alaska in 1996. J. Volcanol. Geotherm. Res. 130, 265–284. doi:10.1016/S0377-0273(03)00302-0

CrossRef Full Text | Google Scholar

Ruppert, N. A., Prejean, S., and Hansen, R. A. (2011). Seismic Swarm Associated with the 2008 Eruption of Kasatochi Volcano, Alaska: Earthquake Locations and Source Parameters. J. Geophys. Res. 116. doi:10.1029/2010JB007435

CrossRef Full Text | Google Scholar

Salah, M. K., Şahin, Ş., and Topatan, U. (2014). Crustal Velocity and Vp/Vs Structures beneath Central Anatolia from Local Seismic Tomography. Arab J. Geosci. 7, 4101–4118. doi:10.1007/s12517-013-1038-7

CrossRef Full Text | Google Scholar

Schmitt, A. K., Danišík, M., Aydar, E., Şen, E., Ulusoy, İ., and Lovera, O. M. (2014). Identifying the Volcanic Eruption Depicted in a Neolithic Painting at Çatalhöyük, Central Anatolia, Turkey. PLoS One 9, e84711. doi:10.1371/journal.pone.0084711

PubMed Abstract | CrossRef Full Text | Google Scholar

Şengör, A. M. C., and Zabcı, C. (2019). The North Anatolian Fault and the north Anatolian Shear Zone. Switzerland: Springer Nature Switzerland AG 2019, 481–494. doi:10.1007/978-3-030-03515-0_27

CrossRef Full Text | Google Scholar

Sparks, S. R. J., Sigurdsson, H., and Wilson, L. (1977). Magma Mixing: a Mechanism for Triggering Acid Explosive Eruptions. Nature 267, 315–318. doi:10.1038/267315a0

CrossRef Full Text | Google Scholar

Stein, R. S. (1999). The Role of Stress Transfer in Earthquake Occurrence. Nature 402, 605–609. doi:10.1038/45144

CrossRef Full Text | Google Scholar

Tank, S. B., and Karaş, M. (2020). Unraveling the Electrical Conductivity Structure to Decipher the Hydrothermal System beneath the Mt. Hasan Composite Volcano and its Vicinity, SW Cappadocia, Turkey. J. Volcano. Geotherm. Res. 405, 107048. doi:10.1016/j.jvolgeores.2020.107048

CrossRef Full Text | Google Scholar

Tezel, T., Shibutani, T., and Kaypak, B. (2013). Crustal Thickness of Turkey Determined by Receiver Function. J. Asian Earth Sci. 75, 36–45. doi:10.1016/j.jseaes.2013.06.016

CrossRef Full Text | Google Scholar

Toda, S., Stein, R. S., Sevilgen, V., and Lin, J. (2011). Coulomb 3.3 Graphic-Rich Deformation and Stress-Change Software for Earthquake, Tectonic, and Volcano Research and Teaching—User Guide. US Geol. Surv. open-file Rep. 1060, 63.

Google Scholar

Toprak, V., and Göncöoḡlu, M. C. (1993). Tectonic Control on the Development of the Neogene-Quaternary Central Anatolian Volcanic Province, Turkey. Geol. J. 28, 357–369. doi:10.1002/gj.3350280314

CrossRef Full Text | Google Scholar

Tosun, H., Zorluer, İ., Orhan, A., Seyrek, E., Savaş, H., and Türköz, M. (2007). Seismic hazard and Total Risk Analyses for Large Dams in Euphrates basin, Turkey. Eng. Geol. 89, 155–170. doi:10.1016/j.enggeo.2006.10.003

CrossRef Full Text | Google Scholar

Ulusoy, İ., Diker, C., Şen, E., Aydın, E., Akkaş, E., Gümüş, E., et al. (2020). Surface Expressions of Morphostructural Features at Hasandağ Stratovolcano on DEM Datasets. Med. Geosc. Rev. 3, 175–191. doi:10.1007/s42990-020-00044-3

CrossRef Full Text | Google Scholar

Walter, T. R., and Amelung, F. (2006). Volcano-earthquake Interaction at Mauna Loa Volcano, Hawaii. J. Geophys. Res. 111. doi:10.1029/2005JB003861

CrossRef Full Text | Google Scholar

Walter, T. R., and Amelung, F. (2007). Volcanic Eruptions Following M ≥ 9 Megathrust Earthquakes: Implications for the Sumatra-Andaman Volcanoes. Geol 35, 539–542. doi:10.1130/G23429A.1

CrossRef Full Text | Google Scholar

Walter, T. R., Wang, R., Zimmer, M., Grosser, H., Lühr, B., and Ratdomopurbo, A. (2007). Volcanic Activity Influenced by Tectonic Earthquakes: Static and Dynamic Stress Triggering at Mt. Merapi. Geophys. Res. Lett. 34. doi:10.1029/2006GL028710

CrossRef Full Text | Google Scholar

Wang, C.-Y., Wong, A., Dreger, D. S., and Manga, M. (2006). Liquefaction Limit during Earthquakes and Underground Explosions: Implications on Ground-Motion Attenuation. Bull. Seismol. Soc. Am. 96, 355–363. doi:10.1785/0120050019

CrossRef Full Text | Google Scholar

Wang, J., Xu, C., Freymueller, J. T., Li, Z., and Shen, W. (2014). Sensitivity of Coulomb Stress Change to the Parameters of the Coulomb Failure Model: A Case Study Using the 2008 M W 7.9 Wenchuan Earthquake. J. Geophys. Res. Solid Earth 119, 3371–3392. doi:10.1002/2012JB009860

CrossRef Full Text | Google Scholar

Wells, D. L., and Coppersmith, K. J. (1994). New Empirical Relationships Among Magnitude, Rupture Length, Rupture Width, Rupture Area, and Surface Displacement. Bull. Seismol. Soc. Am. 84, 974–1002.

Google Scholar

Wilson, L. (1980). Relationships between Pressure, Volatile Content and Ejecta Velocity in Three Types of Volcanic Explosion. J. Volcanol. Geotherm. Res. 8, 297–313. doi:10.1016/0377-0273(80)90110-9

CrossRef Full Text | Google Scholar

Yalcin, A., Gokceoglu, C., and Sönmez, H. (2008). Liquefaction Severity Map for Aksaray City Center (Central Anatolia, Turkey). Nat. Hazards Earth Syst. Sci. 8, 641–649. doi:10.5194/nhess-8-641-2008

CrossRef Full Text | Google Scholar

Young, N. K., and Gottsmann, J. (2015). Shallow Crustal Mechanics from Volumetric Strain Data: Insights from Soufrière Hills Volcano, Montserrat. J. Geophys. Res. Solid Earth 120, 1559–1571. doi:10.1002/2014JB011551

CrossRef Full Text | Google Scholar