A series of Miocene-to-recent calderas, stratovolcanoes and monogenetic fields that constitute the Central Anatolian Volcanic Province (CAVP) (Innocenti et al., 1982) are situated within the Central Anatolian Fault Zone (CAFZ) (Koçyiğit and Beyhan, 1998) (Figure 1). The CAFZ represents a 730 km intracontinental transform fault zone, with ensuing ENE-WSW extension resulting in the regional intra-plate volcanism of the CAVP (Koçyiğit and Beyhan, 1998; Sen et al., 2004). Within the CAFZ, the Erciyes pull-apart basin (EPB) is formed by a S-shaped horst and graben structure creating an extensional depression in which Erciyes Dağ is located (Koçyiğit and Beyhan, 1998).

Erciyes Dağ is an active Pliocene-Quaternary age stratovolcano, with its last known eruption in 6880 BC (Kürkcüoglu et al., 2004). With a summit height of 3,917 m a.s.l., it has a prominence over the surrounding topographic plateau of ∼2,800 m (Gazioğlu et al., 2004). Located ca 25 km south of Kayseri city centre and in the centre of a 14 km by 18 km wide Pleistocene caldera, the volcano towers above the Metropolitan Municipality of Kayseri. With a population of ∼1 million, the southern municipal boundaries reach within a few km of the volcano’s edifice. The Dundarli-Erciyes fault segment of the larger Ecemis fault zone intersects Erciyes’ edifice and the locations of the most recently formed domes along the fault segment are suggestive of magma-tectonic interactions (Koçyiğit and Erol, 2001; Sen et al., 2004). Koçyiğit and Erol (2001) proposed the coeval formation of the EPB alongside Erciyes Dağ’s edifice, with associated crustal thinning encouraging the migration of magma. Seismic data of the crustal architecture of Central Anatolia indicate an overall increase in compressive and shear wave velocities with depth (Salah et al., 2014).

Erciyes Dağ evolved over two principal stages: the initial andesitic-basaltic Koç Dağ stage, and the most recent andesitic-dacitic Erciyes stage (Şen et al., 2003). At the end of the Erciyes stage rhyodacitic domes were extruded along the Dundarli-Erciyes fault segment (Dikkartin Dağ and Perikatin Dağ). Wider geodynamic studies of magmatic sources within Central Anatolia advocate a connection to the ongoing collision of the Eurasian and Arabian tectonic plates (Pasquare et al., 1988) and trace element contributions in magmas from Erciyes Dağ are ascribed to lithospheric assimilation from past subductive periods (Kürkçüoglu et al., 1998). Petrological studies (Dogan et al., 2011, 2013) indicate the presence of two interacting magma reservoirs which fuel eruptions at the volcano: 1) an upper-crustal reservoir at depths of 4–10 km and 2) a deeper-seated (>15 km depth) reservoir that feeds magma to the shallower reservoir (Dogan et al., 2011). The currently preferred petrogenic model for Erciyes describes the flux of parental basaltic melts into the mid to lower crust, at which point they evolve in composition, enriching in silica via processes of crustal assimilation and differentiation (Dogan et al., 2013). These more evolved melts then rise to shallower crustal depths of 4–10 km, at which point they mix with existing entrained crystals and melts within the shallow reservoir (Dogan et al., 2013). This petrogenic model concurs with the concept of transcrustal magma systems occupied by crystal mush and mobile magma (Cashman et al., 2017).

Subsurface mass redistribution and deformation modify the crustal density distribution Δρ(x, y, z). As a result the gravitational potential ϕ _g changes and we quantify this in our models by solving the following differential equation (Cai and Wang, 2005): ∇ 2 ϕ g = − 4 π G Δ ρ x , y , z (1) where G is the universal gravitational constant. The problem is mathematically closed by imposing Dirichlet boundary conditions of zero at infinity.

Here we are concerned with the variations in the vertical component of gravity that is typically measured in volcano gravimetric surveying such that: Δ g x , y , z = − δ ϕ g δ z . (2)

Three source terms (Figure 2) quantify the different contributions to the redistribution of subsurface density (Bonafede and Mazzanti, 1998; Zhang et al., 2004; Currenti et al., 2007): Δ ρ x , y , z = − u ⋅ ∇ ρ 0 + Δ ρ i − ρ 0 ∇ ⋅ u . (3) where u represents the displacement field, ρ ₀ denotes the density of the embedding medium, and Δρ _i is the density change from magma injection into the source.

The first term (denoted Δg₁ hereafter) arises from the shifting of density boundaries within the crust. The second term (Δg₂) quantifies density variations due to the change in mass in the source volume and is composed of two components: 1) the density change resulting from the input of new mass into the reservoir at constant volume V and the compression of resident magma and 2) the density change accompanying the volume expansion ΔV of the reservoir. The first component is dependent on reservoir compressibility β _r = β _c +1/ρ × Δρ/ΔP, where β _c is the compressibility of the encasing host rock, while the second component is dependent on the density change induced by the deforming reservoir and the replacement of mass surrounding the reservoir. The third term (Δg₃) reflects modulation of the changes in source volume by host-rock compressibility (Bonafede and Mazzanti, 1998).

To obtain residual changes in the vertical component of gravity (Δg _r ), one needs to also account for the “Free-Air” gravity change Δg₀ that accompanies vertical surface deformation: Δ g 0 = γ w (4) where γ is the Free-Air gradient (−308.6 µGal m⁻¹; 1 µGal = 10^–8 m s⁻²) and w is the vertical change in the position of the observation point.

All four terms contribute to the residual gravity variation: Δ g r = Δ g 0 + Δ g 1 + Δ g 2 + Δ g 3 (5)

In this study we only consider elastic deformation in the form of Hooke’s law which relates stress σ and strain ϵ via: σ = 2 μ ϵ + λ t r a c e ϵ I (6) ϵ = 0.5 ∇ u + ∇ u T (7) where μ and γ are the Lamé parameters and I is the identity matrix.

Eastward, northward and vertical displacement vectors u, v and w, respectively, are derived directly from the model from which total displacements are calculated.

Ground displacement and gravity change data must be modelled jointly and simultaneously to quantify the different contributions to changes in the gravitational potential from subsurface processes.

We solve the above equations numerically using Finite Element Analysis in COMSOL Multiphysics v5.4 for a suite of 2D axisymmetric and 3D forward models to study the effect of crustal mechanical heterogeneity and topography on measured ground displacements and gravity changes from rejuvenation of Erciyes Dağ’s magmatic reservoir. Model parameters and symbols utilised throughout this study may be found in Table 1.

Sets of different 2D axisymmetric and 3D models are developed to systematically explore the influence of crustal heterogeneity and topography on source parameter changes (generic models) and minimum detectable mass and volume fluxes (minimum detectability models). Models with a flat free surface are labelled “FLAT,” while models implementing a topography are labelled “TOPO”: 1) Base models (BM) including FLAT models (FM) and Layered FLAT models (LFM), and 2) TOPO models (TM) including layered TOPO models (LTM) and non-layered TOPO models (NLTM). 2D models are covered in the Supplementary Material. Here we focus on the development and analysis of 3D models with model characteristics summarised in Table 4.

We determine mass and volume changes upon reservoir recharge for given source pressurisation (ΔP) and source density change (Δρ) at the detectability limits of joint gravimetric and geodetic field surveys. We set the limits to 0.01 m for vertical surface displacements and 0 ± 5 μGal (1μ = 10^–8 m s⁻²) for residual gravity changes based on the typical sensitivity of field instrumentation and survey protocols as well as constraints from remote sensing of volcanoes in the CAVP (Biggs et al., 2021). The minimum detectable volume flux into the reservoir is given by McTigue (1987): Δ V = α π μ Δ P (12) where α is the source radius and μ is the shear modulus. The associated reservoir minimum mass change is calculated from Δ M = Δ V ρ m (13) using the relationship between mass change, density and volume change.

All results from generic models are reported in Table 4. Here we focus on presenting the differences in results by comparing different modelling approaches as a function of increasing complexity.

Value pairs of ΔP and Δρ that produce mass and volume changes upon magma reservoir recharge at the minimum detectability limit of surface uplift (0.01 m) and residual gravity changes (±5 µGal) are reported in Table 5. To illustrate, for the assumptions underpining the LTM we derive a pressure change of 0.25 MPa and concurrent density change of ±0.067 kg m⁻³ at the detectability limit.

The LFM predicts the smallest ΔP values with the TM predicting the highest. Δρ values vary between ±0.067 kg m⁻³ in the LFM and ±1 kg m⁻³ in the TM. The LFM and TM bracket the range in values at the lower and upper end, respectively. The deduced reservoir mass and volume changes vary accordingly across the different models. The NLTM predicts >50% larger mass and volume changes at the detectability compared to the LTM. Results for the “simplest” model (i.e., a FM) fall between the values predicted for the most complex models (i.e., LTM and NLTM). Figure 7 shows the contrast in resolvable source changes in the TM, LTM and NLTM with respect to those derived from the FM.

Our study shows the limited applicability of traditional homogeneous elastic half-space models to inform processes at Erciyes Dağ. Though usually preferred in volcano geodetic modelling, the elastic homogeneous half-space Mogi model (Mogi, 1958) only provides accurate solutions for source depth: source radius ratios of >5. At Erciyes Dağ a maximum source radius for Erciyes Dağ of 1.4 km would be permissible for a source depth of 7 km derived from petrological data (Dogan et al., 2013). Consequently, the 3.1 km source radius calculated from previous eruption volumes is more than twice larger than the maximum permissible value, causing the Mogi solution to underpredict both vertical and horizontal deformation at the surface for given source pressurisation. This agrees with our original hypotheses and matches the results of previous studies (e.g., Battaglia and Segall, 2004), providing support to prioritise numerical models over analytical models to study ground deformation and gravity variations at Erciyes Dağ under the assumption that the reservoir source depth: source radius ratio is ≪5.

Subsurface heterogeneities fundamentally modulate the stress vs. strain relationship, affecting surface displacements (Gudmundsson and Brenner, 2004; Geyer and Gottsmann, 2010; Gottsmann et al., 2020) and gravity variations by either muting or amplifying strain depending on the subsurface distribution of mechanically stiff and weak lithologies. The impact of subsurface heterogeneities on changes in gravity is that they modulate the subsurface density distribution expressed by the u and ∇u terms of Eq. 2. This effect is demonstrated by the >100% increase in w, Δg₁ and Δg₃ components and the >700% increase in Δg₂ in the LFM compared to the FM. The extraordinary change in the magnitude of Δg₂ between homogeneous and heterogeneous models is attributed to the more pronounced expansion of a reservoir embedded in mechanically weaker rocks and the associated shifting of density boundaries permitting the influx of additional mass into the reservoir volume. Most notably, the resultant residual gravity change is 14 times greater in the LFM compared to the FM.

In our models, the combination of a mechanically weak volcanic edifice and upper-crustal mechanical weaknesses exert the greatest modulation of surface strain and the gravity field from upper-crustal reservoir pressurisation (Figures 5, 6). It thus follows that the magnitudes of the Δg₁ and Δg₃ terms correlate with higher degrees of mechanical complexities; e.g., the number of abrupt changes in the subsurface density distribution.

In comparing predicted residual gravity changes from FM with those from LTM and NLTM for the set of generic model parameterisations (Table 3), their magnitude is greatly increased primarily by amplified Δg₂ contributions. Currenti et al. (2007) analysed the effects of subsurface heterogeneities using 2D axisymmetric models for a smaller source radius (1 km) and larger pressure changes (100 MPa). Whilst their observations are analogous to ours in relation to the amplification of gravity changes and displacements, their calculated magnitudes are smaller owing to a diminished Δg₂ term due to the mass addition to a smaller source with a greater overburden producing smaller crustal strain.

The displacement of density boundaries between the source and host rock, along with the free-air effect, permit significant reservoir mass changes at net source density changes of ∼0 kg m⁻³. This signifies the importance of considering contributions to gravity changes from a heterogeneous subsurface in volcanic areas. Modelling approaches relying on the assumption of a mechanically homogeneous subsurface would need to attribute residual gravity changes almost exclusively to source density changes with major implications for the interpretation of causative subsurface processes amid volcano uplift. Judging from the Δρ values reported in Table 5 for a FM compared to a LFM, around one order of magnitude larger reservoir density changes would be deduced under the assumption of medium homogeneity.

Implementing a gradual change in mechanical properties instead of discrete changes in the form of layers has a number of implications for modelling approaches and the interpretation of results. First, non-layered models solve for <50% of the degrees of freedom of layered models at a much reduced computational cost. This cost efficiency could be compensated, for example, by the inclusion of a greater number of complexities, such as a higher resolution topographic surface. Second, in our study generic non-layered models produce vertical and total displacements that are ∼50% of the magnitude of layered models due to the amplification of surface strain by low rigidity discrete layers in the upper crust (Geyer and Gottsmann, 2010; Ronchin et al., 2015). Third, whilst reducing computational cost, non-layered models predict significantly different magnitudes of gravity changes compared to layered models. The large diminution of the Δg₂ contribution in non-layered models results in a residual gravity change over 5 times smaller than in layered models. Because the crust is neither composed solely of discrete density layers, nor is it best described by a single density gradient (Zhu et al., 2019), the layered and non-layered models represent end-member scenarios of the crustal architecture. Source property changes deduced from the evaluation of gravity changes and surface displacements using either mechanically heterogeneous layered or non-layered models likely map the uncertainties behind data interpretations.

The inclusion of topography in models to study subsurface processes at high-prominence volcanoes results in decreased magnitudes of surface displacements compared to models ignoring topography for given source changes. In general we find that edifice prominence and magnitudes of ground displacement are negatively correlated. This is in agreement with previous studies (e.g., Cayol and Cornet, 1998; Williams and Wadge, 1998; Ronchin et al., 2015), and is due to the addition of a crustal mechanical domain that influences the subsurface stress and strain relationship. The effect of displacement modulation decreases with increasing source depth, particularly with respect to vertical displacement, and hence has a major influence on gravity source terms 1 and 3 given the 53% reduction in w in a FM compared to a TM (Table 4).

Muting the effect of crustal straining by source pressurisation, coupled with the mass-related gravitational attraction being chiefly determined by the distance between the computational surface and the source (Charco et al., 2009), volcano prominence has a major influence on deduced source mass and volume changes from monitoring observables. Whilst subsurface heterogeneities influence the magnitude of modeled anomalies, topography modifies both the magnitude and the spatial pattern of surface displacements and gravity changes. Prominence contributes to a flattening of the peak magnitude of both observables, with the largest impact at Erciyes Dağ noted directly over the summit where the slope angle is greatest (Figures 5, 6). These findings broadly corroborate results from axisymmetric models presented in Cayol and Cornet (1998) and Charco et al. (2007) and attest to the importance of considering topography for interpreting data from large-prominence volcanoes.

A significant mass addition of between 10⁹ and 10¹⁰ kg to the reservoir feeding Erciyes Dağ is possible without detection at the surface. The range of deduced values is model dependent with mechanically homogeneous models requiring larger source changes than layered models at the detectability limits. These disparities are explained by the inclusion of subsurface layers which result in an increase in gravity changes and surface displacements for given source changes (Currenti et al., 2007). This amplification reduces the minimum detectable reservoir mass change as a smaller mass input is required to achieve surface changes at the detection limit.

We have used available multi-parametric data for Erciyes Dağ to constrain reservoir shape and location as well as crustal mechanical heterogeneity to the best of our knowledge. Implementing a 1D heterogeneity profile based on low-resolution seismic tomography data of Central Anatolia augmented by shallow crustal data from Mt. St. Helens and Soufrière Hills Volcano likely affect the accuracy of the results. However, with broadly similar near-surface geologies, upper-crustal and edifice mechanical characteristics of the three volcanoes, we do not consider our results to be fundamentally flawed. A much greater bias of model results are expected from the assumption of source sphericity and crustal elasticity in response to magma rejuvenation. Thermomechanical effects (Hickey et al., 2015; Morales Rivera et al., 2019; Gottsmann et al., 2020) are expected to affect results with the implication that inelastic accommodation of magma would render the derived mass and volume changes lower-bound estimates. The reservoir dimensions are such that they satisfy both available petrological, volumetric and volcanological constraints, rendering a spherical geometry as the most plausible. An alternative geometric representation of the reservoir is a vertically extensive spheroid (e.g., a prolate ellipsoid). Such a reservoir geometry yields smaller vertical displacements, resulting in a more nuanced free-air effect and smaller magnitudes of Δg₁ and Δg₃ contributions. As a result, larger source volume changes and larger mass changes compared to those derived for a spherical reservoir are required to be detectable. However, a reservoir of prolate ellipsoidal geometry matching the volume estimate of 124 km³ would need to extend to depths well beyond those constrained by petrological evidence. Hence, a prolate geometry would also imply a smaller reservoir volume. To simultaneously address the potential overestimation of reservoir volume and a non-spherical reservoir geometry, we developed a NLTM for a prolate reservoir of major-semi axis = 3.1 km and semi-minor axes of 2.2 km (matching a reservoir volume of ∼62 km³, i.e., ∼50% of the derived volume in Eq. 1). This model requires a ∼50% larger source volume change and a ∼45% higher source density change at the detectability limit compared to the NLTM for a spherical source of 124 km³. This indicates that, given the data currently available for Erciyes Dağ, model results are less sensitive to reservoir geometry.

Our models ignore potential volcano-tectonic linkages (Figure 9). The Dundarli-Erciyes fault segment which intersects Erciyes’ edifice and caldera boundary faults may 1) encourage the migration of magma from the source to the surface along the fault systems and 2) permit volcano-tectonic stress interactions. The location of the most recently formed domes along the fault segment and its intersections with the caldera boundary fault (Şen et al., 2003) matches observations of faults providing pathways for fluid and magma migration (Miller et al., 2017). In all our models, the volcanic reservoir is located directly beneath the edifice. Whilst informed by petrological data, the lack of geophysical data preclude further constraining of the location of the reservoir. A reservoir positioned further from the edifice would be expected to require smaller mass and volume changes at the detectability limit due to the reduced influence of topography. To test this we developed a NLTM with the reservoir’s centroid located 5 km south of the edifice as shown in Figure 9. The predicted residual gravity changes and total displacements for generic source change parameters reported in Table 3 were up to 124 and 117% larger, respectively, compared to the results from the centroid’s position directly beneath the edifice. This implies mass and volume changes at the detectability limit that are a factor of ∼2 smaller than derived from models with a reservoir centered beneath the summit.