Minimum Detectable Mass and Volume Fluxes During Magmatic Recharge at High Prominence Volcanoes: An Application to Erciyes Dağ Volcano (Turkey)

Magma reservoir recharge is widely recognised as a precursor of eruptive activity. However, the causative relationships between reservoir rejuvenation and surface observables such as gravitational potential field changes and ground deformation are still poorly understood. At intermediate and silicic intra-plate volcanoes where crustal mechanical heterogeneity combined with high-prominence are expected to fundamentally affect the crustal stress and strain relationship, protracted period of repose and absence of monitoring data raise questions about the detectability of magma recharge. Here we report results from integrated geodetic forward modelling of ground displacements and gravity changes from reservoir recharge at Erciyes Dağ, a large prominence (∼2,800 m), yet poorly studied, stratovolcano of the Central Anatolian Volcanic Province in Turkey. The most recent eruption at ∼7000 BC, close proximity to the Kayseri Metropolitan Area and absence of dedicated volcano monitoring set a precedent to explore stealth magmatic processes at the volcano. Using finite element analysis we systematically explore the influence of subsurface mechanical heterogeneities and topography on surface deformation and gravity changes from magmatic recharge of Erciyes Dağ’s reservoir. We show that whilst crustal heterogeneity amplifies ground displacements and gravity variations, the volcano’s substantial prominence has the opposite effect. For generic reservoir pressure and density changes of 10 MPa and 10 kg m−3 predicted vertical displacements vary by a factor of 5 while residual gravity changes vary by a factor of 12 between models ignoring topography or mechanical heterogeneity and those that do not. We deduce reservoir volume and mass changes of order 10–3 km3 and 1010 kg, respectively, at the detectability limit of conventional surveying techniques at the volcano. Though dependent on model assumptions, all results indicate that magma recharge at Erciyes Dağ may go undetected at fluxes 1) sufficient to maintain an active reservoir containing eruptable magma and 2) similar to those reported for intermediate/silicic volcanoes with repose times of 100–1,000s of years (e.g., Parinacota) and persistently active mafic volcanoes such as Mt. Etna and Stromboli. Our findings may be utilised to inform integrated geodetic and gravimetric monitoring at Erciyes Dağ and other large prominence silicic volcanoes and could provide early insights into reservoir rejuvenation with implications for the development of disaster risk reduction initiatives.


INTRODUCTION
Joint ground deformation and microgravity observations are important methods for assessing subsurface processes at active volcanoes including during unrest and the lead-up to eruptive activity (Rymer and Brown, 1989;Carbone et al., 2007;Battaglia et al., 2008;Currenti, 2014;Carbone et al., 2017;Fernández et al., 2017;Gottsmann et al., 2020). Coupled with data modelling using both analytical and numerical approaches, these techniques allow the delineation of spatio-temporal variations in density and mass in sub-volcanic plumbing systems. Most analytical geodetic models rely on the concept of a causative source embedded within an isotropic, homogeneous, elastic half-space beneath a flat free surface (Mogi, 1958;Hagiwara, 1977). However, it has been acknowledged that such simplifications may lead to biased results on source parameters because volcanic regions have 1) substantial topographical relief, with some volcanic edifices showing significant >2,000 m topographic prominence (Cayol and Cornet, 1998;Williams and Wadge, 1998) and 2) are characterised by complex crustal heterogeneities due to the nature of geodynamic and eruptive processes underpinning volcanism (Gudmundsson and Brenner, 2004). Currenti et al. (2007) highlight a negative correlation between topographic prominence and detectable ground displacements and residual gravity variations. Medium heterogeneity fundamentally affects the subsurface stress vs. strain relationship and thus the resultant deformation measured at the surface (Bonaccorso et al., 2005;Hickey et al., 2016). The magnitude of surface displacement and thus gravity variations may be amplified, or muted, by mechanically soft and stiff lithologies, respectively, with amplification effects intensified at shallower depths where crustal layers are typically softer (Geyer and Gottsmann, 2010). Topographic and mechanical complexities necessitate the use of numerical approaches which come at substantial computational costs and are hence often ignored in favour of analytical approaches. This stands in contrast to observations, particularly at steep-sided stratovolcanoes, which show a significant number of eruptions that occur with no prior surface deformation (Biggs et al., 2014). Here we focus on Erciyes Dağ in Central Anatolia, Turkey, an active and large prominence volcano that lacks dedicated volcano monitoring. Hence little to nothing is known about subsurface processes currently operating at the volcano and their relation with surface observables. Exploring both analytical and numerical modelling, this study aims to determine reservoir mass and volume changes at Erciyes Dağ from minimum resolvable surface deformation and gravity changes. The hypothesis underpinning this work is that magma reservoir rejuvenation at Erciyes Dağ may occur at below the detectability limit of spatiotemporal geodetic and gravimetric anomalies.

Central Anatolian Volcanic Province
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ğ
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(Dogan et al., , 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).  Dogan et al. (2013). The Anatolian plate is bounded by the Northern and Eastern Fault Zones between the colliding Eurasian and Arabian plates. (B) Regional volcano-tectonic map of Erciyes Dağ after Koçyiğit and Erol (2001). The volcano is located 25 km south of the city centre of Kayseri and within the Erciyes pull-apart basin (EPB

Potential Field and Ground Deformation Analysis
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): 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: 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): where u represents the displacement field, ρ 0 denotes the density of the embedding medium, and Δρ i is the density change from magma injection into the source. The first term (denoted Δg 1 hereafter) arises from the shifting of density boundaries within the crust. The second term (Δg 2 ) 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 3 ) 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 0 that accompanies vertical surface deformation: where c is the Free-Air gradient (−308.6 µGal m −1 ; 1 µGal 10 -8 m s −2 ) and w is the vertical change in the position of the observation point. All four terms contribute to the residual gravity variation: In this study we only consider elastic deformation in the form of Hooke's law which relates stress σ and strain ϵ via: where μ and c are the Lamé parameters and I is the identity matrix.
FIGURE 2 | Schematics of the principal contributions to gravity changes measured at volcanoes from changes in the subsurface density distribution produced by: (A) the displacement of density boundaries (including the free surface, the source walls, and crustal density boundaries), (B) the source mass change, and (C) the volume change of the reservoir itself modulated by crustal compressibility (after Bonafede and Mazzanti (1998) and Currenti et al. (2007)).
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 750063 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.

Model Development and Parameterisation
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.

Subsurface Mechanics
Geophysical data on the subsurface architecture of Erciyes Dağ are largely absent. Seismic tomography data by Salah et al. (2014) provide low resolution 1D p-wave velocities (V p ) and p-wave/swave velocity ratios (V p /V s ) across Central Anatolia for depths of 4, 12, 25 and 40 km which we use to parameterise subsurface mechanical heterogeneity ( Figure 3). Owing to the absence of seismic data for the uppermost (<4 km depth) Central Anatolian crust, we use seismic data from the Cascades around Mt. St. Helens (Kiser et al., 2016) as a geologically plausible alternative due to similar seismic velocity distribution at depth >4 km (e.g., at 4 km depth V p and V p /V s data in Central Anatolia are only <1.5% larger than in the Cascades). We use the seismic data to derive elastic properties following empirical equations presented in Brocher (2005): Depth dependent Poisson's ratios (]) were determined from V p /V s ratios: to derive the dynamic Young's modulus E d from: We convert E d to static Young's moduli E s for appropriate pressures and temperatures following Eissa and Kazi (1988) such that E s 0.5E d . Parameters are computed for each depth segment of the crust shown in Figure 3, for 0 > z > -34 km (the depth of the Moho), with values reported in Table 2. Note that in all models, depth is expressed by negative values of z. For clarity, in the text we refer to depth by positive values. Median average values of subsurface mechanics used in generic models are shown in Table 3. Elastic properties for Erciyes Dağ's edifice were taken from Montserrat (Young and Gottsmann, 2015) due to broadly similar eruptive products of dome-building andesite volcanism at both volcanoes. In all derivations, effects of temperature variations in the subsurface on elastic properties are neglected.

Magma Reservoir
In the absence of any geophysical constraints, we use petrological data in combination with calculations of the magma chamber volume to estimate source parameters such as depth and geometry. Dogan et al. (2013) used amphibole geobarometry to define an upper crustal magma reservoir located at 4-10 km depth beneath Erciyes Dağ within a homogeneous crust. We recalculated the barometric data to derive reservoir depth account for the subsurface density distribution. This provided a centroid depth of the reservoir of z −7 km beneath the surface. The total volume of a magma reservoir (V m ) with magma compressibility (β m ), located in country rock with in-situ tensile strength (T 0 ) and average compressibility (β c ), is related to the volume of ejected material (V e ) for a given eruptive phase (Browning et al., 2015) via: FIGURE 3 | Schematic of subsurface material properties (crustal density ρ c and crustal static Young's modulus E s ) of the layered TOPO model (LTM) derived from 1D seismic tomography data of central Anatolia (Salah et al., 2014) using Eqs 8-10. In the absence of shallow crustal data for the volcano and surrounding areas, we use seismic data from the Cascades around Mt St Helens (Kiser et al., 2016). Edifice properties are taken from Soufrière Hills Volcano (Young and Gottsmann, 2015). An idealised topographic relief of Erciyes Dağ is shown for simplicity rising ∼3 km from the surrounding plateau (z 0). Digital elevation data is used to parameterise the topographic surface in the Finite Element models. The magma reservoir is represented by a sphere with its centroid located directly beneath the summit of the volcano at a depth of z −7 km.
Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 750063 We calculate an average value of β c 3(1−2]) E s 1.78 × 10 -11 Pa −1 for the reservoir depth range and set β m 7.6 × 10 -10 Pa −1 to account for the explosive nature of Erciyes' dome-forming eruptions (Voight et al., 2010). A dense rock equivalent (DRE) V e 0.958 km 3 for Erciyes' most recent dome-forming eruption at Dikkartin Dağ (Şen et al., 2002) gives V m 123.5 km 3 . In the absence of any information that would indicate otherwise, a spherical reservoir of radius α 3.1 km is the most plausible representation of the reservoir given the volumetric and petrological constraints. In all models the surface projection of the centroid is at 38°31′54″N and 35°26′48″E.

Model Suite
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.

3D Base Models
An initial 3D FM is created to benchmark against a 2D axisymmetric FM. An 80 km * 80 km * 68 km computational domain is created in order to prevent artefacts within the numerical solution inflicted by boundary effects. The base of the model is set at z −34 km, representing the Moho depth, with the free surface, representing the volcanic plateau, set at z 0 km. Median average values for mechanical properties reported in Table 3 are applied to invoke subsurface mechanics. Dirichlet boundary conditions are set to zero on the outer faces, a fixed base is added, and roller conditions are applied to the lateral faces of the domain for z < 0 km. Coefficient form Partial Differential Equations (PDE) representing the Δg 1 , Δg 2 and Δg 3 terms are solved for the edifice, crust and source to calculate gravitational potential variations at the surface. The Δg 0 term is accounted for during post-processing. A 3D LFM is then compared to a 3D FM to resolve the effect of subsurface density variations within the crust. The LFM incorporates eight layers with material properties defined in Table 2.

3D Layered TOPO Models
Building on the LFM we implement a 3D topography using 90 m SRTM data to generate a digital elevation model (DEM) of the surface to obtain an LTM. The DEM is imported and introduced as a parametric surface with a relative tolerance of 0.001 and a maximum number of 150 knots to optimise resolution of the surface whilst precluding excessive computational cost. The source is implemented at the coordinates as in the BM and sits directly beneath the edifice at a centroid depth of 7 km beneath the plateau. Parametric surfaces representing Erciyes' subsurface layers are employed, along with a cylinder of radius 10 km extending to a depth of 1 km around the volcanic edifice for meshing purposes (Figure 4). The PDEs are solved to each subsurface layer individually. Meshing follows a structured approach whereby mesh refinement is applied to the cylindrical domain encompassing the edifice and the source (Hickey et al., 2015). Through a computational cost-benefit analysis, it was determined that higher resolution, low-cost models may be obtained by adding a point to the topographic surface at the centre of the edifice and creating a custom mesh around it. Maximum and minimum element sizes of 20 and 5 m, respectively, are defined with a maximum element growth rate of 1.1 and a curvature factor of 0.4. The effect of topography in integrated geodetic modelling is evaluated by comparing results from the LTM against those from the LFM. Additionally, the LTM results are contrasted with those from the TM to assess the effect of subsurface layering.

3D Non-Layered TOPO Models
Whilst the mesh refinement technique described above prevents high computational cost, the expense is still large. To further reduce cost, 3D NLTMs are produced, decreasing the required quantity of mesh elements due to fewer subsurface layers. Polynomial fits to ρ c , E s , and ] data presented in Table 2 are TABLE 2 | Mechanical properties assigned to subdomains in layered models based on seismic velocities (V p V s ). In the absence of local seismic data for the edifice and uppermost portions of the crust, data from Soufrière Hills volcano Young and Gottsmann (2015) are used for z > 1 km and from the Cascades Kiser et al. (2016) for 1 ≥ z ≥ −3 km. Values for z < -3 km are derived from data presented in Salah et al. (2014). The plateau surrounding Erciyes Dağ is at z 0 km in the model, thus the edifice extends to a height of 2,917 m. See Eq. 8 and Table 1 for further information on the derivation of the values.  used to parameterise the NLTM (Figure 2. The model domain geometry is set up as before minus the parametric surfaces defining the layers (Figure 4). The cylindrical domain of radius 10 km remains around Erciyes' edifice due to the greater disparity between the edifice and subsurface crustal properties, as do the material properties assigned to the edifice and source. Model physics and a structured mesh as described above are added to the domain, and Δg 1 , Δg 2 and Δg 3 terms are solved for the edifice, crust and source. Results are compared to those of the LTM to determine the effect of discrete vs. gradual subsurface layering on gravity variations and surface deformation. The NLTMs provide insight on the gravity change contributions from shifting density boundaries in a layered versus non-layered crust. Moreover, depending on acceptable uncertainties within the results, the penalty in terms of accuracy between models may be determined given that the crust is neither fully density layered, nor does subsurface density gradually vary throughout. The LTM solves for ∼5.5 million degrees of freedom, whereas the NLTM solves for ∼2.3 million degrees of freedom and solves ∼50% faster.

Mass and Volume Fluxes at the Detectability Limit
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 −2 ) 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): where α is the source radius and μ is the shear modulus. The associated reservoir minimum mass change is calculated from using the relationship between mass change, density and volume change.

Generic Models
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. FIGURE 4 | Schematics of 3D model setups for the heterogeneous layered and non-layered TOPO models (LTM and NLTM) including applied domain boundary conditions. Whilst the subsurface of the layered model contains discrete mechanical layers throughout the crust, the nonlayered model comprises gradually varying mechanical properties. The two models denote end-member implementations of a heterogeneous mechanical crustal structure beneath Erciyes Dağ.

TM Versus LTM
Subsurface mechanical layering increases magnitudes of all model outputs compared to those from implementing homogeneous mechanics. Total displacement and w increase 2.8 fold between the TM and the LTM, with u and v increasing 2.9 and 3 times, respectively ( Figure 5). The maximum magnitude of Δg 1 increases by a factor of 2.7 in the LTM whereas the minimum Δg 1 value increases 2.3 fold. The Δg 2 term is increased 8.6 times and the Δg 3 contribution is increased by a factor of 3.3 in the LTM compared to the TM. The resulting Δg r value increases 13.3 fold in the LTM ( Figure 6).

LTM Versus NLTM
Approximating subsurface mechanical heterogeneity by gradual changes in properties mutes the magnitudes of all results compared to models implementing discrete changes in subsurface mechanics. Maximum magnitudes of u and v are decreased by 39 and 38%, respectively, in the NLTM, with total and w displacements decreased by 48% between the LTM and NLTM ( Figure 5). The NLTM predicts maximum values of Δg 1 that are 36% smaller than the LTM and minimum values that are 54% smaller. The Δg 2 contribution is reduced by 74% in the NLTM, with the Δg 3 contribution reduced by 41%. These changes result in a reduction of the Δg r values by 81% in the NLTM compared to the LTM ( Figure 6).

Minimum Detectability Models
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 −3 at the detectability limit. The LFM predicts the smallest ΔP values with the TM predicting the highest. Δρ values vary between ±0.067 kg m −3 in the LFM and ±1 kg m −3 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.

DISCUSSION
At long-quiescent volcanoes with no monitoring data, the analysis and interpretation of synthetic results produced by models such as those developed in this study is challenging. Results are highly sensitive to assumptions underpinning modelling approaches (Tables 4, 5). To better understand the impact of topography and mechanical heterogeneity on deduced mass flux and volume changes due to magma recharge at Erciyes Dağ, we disentangle the individual contributions to surface gravity variations and displacements from our models.

Application of Analytical Models
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.

The Effect of Subsurface Heterogeneities
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 1 and Δg 3 components and the >700% increase in Δg 2 in the LFM compared to the FM. The extraordinary change in the magnitude of Δg 2 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 1 and Δg 3 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 2 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 2 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 −3 . 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 nonlayered 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 2 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 Effect of Topography
The inclusion of topography in models to study subsurface processes at high-prominence volcanoes results in decreased magnitudes of surface displacements compared to models FIGURE 7 | Percentage difference in reservoir parameter changes from TM, LTM and NLTM with respect to the FM to produce minimum detectable mass and volume changes upon magma recharge (see also 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.

Minimum Detectable Mass Intrusions and Volume Changes at Erciyes Dağ
A significant mass addition of between 10 9 and 10 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 . This amplification reduces the minimum detectable reservoir mass change as a smaller mass input is required to achieve surface changes at the detection limit.

Simple Vs. Complex Models
Whilst the minimum detectable ΔP and Δρ values in models employing both a flat surface and a homogeneous subsurface (i.e., a FM) are distinct from the results of more complex models, associated mass and volume changes for the FM fall between results for the most complex models (LTM and NLTM). One could be tempted to argue that the simplest model yields results that do not differ too much from the ones obtained by the more complex models (LTM and NLTM). In our study the simplest model (FM) already implements Median subsurface mechanical parameters derived from seismic studies, and does not utilise generic values that are often applied in analytical geodetic modelling. We would expect differences in results between more complex models and simpler models to become larger with increasing simplifications of the "true" subsurface mechanics. Source parameters from inverse models implementing a 3D mechanical architecture are significantly different to those from isotropic, homogeneous and elastic half-space models (Hickey et al., 2016). The marked differences between results from FMs compared to more complex models become evident for plausible (yet generic) changes in source parameters for surface observables significantly above the detectability thresholds (Figures 5, 6) where predicted ground upift varies by a factor of ∼5 and residual gravity changes differ by more than one order of magnitude for the same source processes.
The principal disparities between the layered and non-layered models are the increases of magnitudes in Δg 1 and Δg 2 terms. Since the Δg 1 term relates to the shifting of subsurface density boundaries, the minimum detectable mass influx upon magma recharge is greater in non-layered models as a result of the removal of discrete boundary layers. It thus follows that the Δg 2 term related to the introduction of new mass increases in non-layered models as a greater mass influx may occur before detection thresholds are reached.

Magma Reservoir Dynamics
The influx of magma at Erciyes Dağ's reservoir is likely modulated by a highly compressible magma mush given the explosive nature of its most recent dome-forming eruptions. Magma reservoirs feeding stratovolcanoes have characteristically high magma compressibilities and distinctively higher rates of eruptions that are not preceded by measurable surface deformation (Biggs et al., 2014) due to a combination of crustal mechanics, local stress fields and high H 2 O-contents (Mastin et al., 2008;Ebmeier et al., 2013).
Notwithstanding the influence of magma compressibility on the detectability limits of ground displacements and gravity changes, the shallow depth and large size of the reservoir modelled in this study create larger displacements and gravity field variations than a similar-sized reservoir at greater depth. This implies that the deduced mass and volume fluxes represent the lower-bound estimates of plausible, yet, undetectable magma rejuvenation.
Minimum magma volume fluxes of between 5 × 10 -3 and 10 -2 km 3 yr −1 are required to sustain large volumes of mobile magma and prevent their solidification (Annen, 2009;Gelman et al., 2013). Whilst these magma fluxes provide a useful reference for the magma fluxes required for caldera-forming eruptions, Erciyes Dağ's most recent eruptions suggest that dome-forming eruptions are more likely to occur in the future. Annen et al. (2015) determined that to prevent complete crystallisation of a volcanic storage region, a time-averaged magma flux exceeding 10 −4 -10 −3 km 3 /yr is required. This emplacement rate would permit the development of a storage region into a potentially active magma reservoir.
These findings relate to the significant source volume increase of 2-5 × 10 -3 km 3 upon magma recharge at Erciyes Dağ without observable surface displacement or residual gravity variations. Under the reasonable assumption that magma reservoir replenishment at Erciyes Dağ occurs below the detection limit of surface uplift velocities of 0.01 m yr −1 (Biggs et al., 2021), the Frontiers in Earth Science | www.frontiersin.org October 2021 | Volume 9 | Article 750063 reservoir may currently sustain eruptable magma if deduced mass changes are taken as annual fluxes. This "stealth" magma flow into the reservoir matches and exceeds reported magma fluxes at several intermediate and silicic stratovolcanoes with repose times of 100-1000s of years (e.g., Parinacota, Soufrière Hills volcano and Mt. Pelée) as well as persistently active mafic volcanoes such as Mt. Etna and Stromboli (Figure 8). This likely attests to fundamentally different geometries and dynamics of magmatic plumbing systems at persistently active mafic volcanoes with high-level magma storage compared to non-mafic volcanoes with mid-to lower crustal magma reservoirs and protracted repose times. It is therefore distinctly possible that Erciyes Dağ has an active magma reservoir fed by magma intrusions that are accommodated by crustal elasticity and high magma compressibility. The presence of eruptable magma may imply that future eruptions may occur at short notice and within a matter of days or weeks (White and McCausland, 2019).

Model Limitations
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 1 and Δg 3 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 3 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 3 , 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 3 . 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  (Branca and Ferrara, 2013), and Stromboli (Allard et al., 1994). The blue dashed line represents the minimum volume flux that is required to retain mobile and eruptable magma within the reservoir (Annen et al., 2015). The purple dashed line represents the minimum volume flux required to produce caldera-forming super-eruptions (Gelman et al., 2013). Minimum detectable magma fluxes for Erciyes Dağ are representative of the two end-member models (LTM and NLTM).
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.

CONCLUSION AND IMPLICATIONS
Reservoir volume and mass changes of orders 10 -3 km 3 and 10 10 kg, respectively, derived for Erciyes Dağ at the detectability limit of integrated geodetic surveys support the hypothesis that magma recharge within a mechanically heterogeneous crust beneath high-prominence volcanoes may go undetected using conventional geodetic and gravimetric monitoring techniques. Ignoring subsurface heterogeneity and topography in models biases the interpretation of gravitational potential field changes and surface displacements. First, Erciyes Dağ's high prominence significantly decreases the magnitude and alters the spatial pattern of surface displacement and gravity variations. Second, mechanically compliant lithologies in the upper crust and edifice amplify ground displacements and gravity changes.
Our findings demonstrate that geodetic and gravity time series data must be jointly modelled and interpreted. This permits the quantification of subsurface mass and volume fluxes upon reservoir recharge which is particularly important at stratovolcanoes where a significant proportion of eruptions occur with no prior ground deformation due to high magma compressibility (Biggs et al., 2014). The deduced amplitude of reservoir replenishment at Erciyes Dağ commensurate with the geodetic detection limit of 0.01 m yr −1 of uplift is potentially sufficient to maintain an active magma reservoir containing eruptable magma. Based on our findings we recommend the implementation of routine joint geodetic and gravimetric monitoring of Erciyes Dağ for disaster risk prevention owing to the volcano's proximity to the densely populated Kayseri Metropolitan Area.

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: Zenodo repository (https://doi.org/10.5281/zenodo.5142185).
FIGURE 9 | Volcano-tectonic map of Erciyes Dağ and surface projections of explored magma reservoirs (to scale). Dikkartin Dağ and Perikartin Dağ, the two most recent eruptive centres of the volcanic complex are located along the Dundarli-Erciyes fault segment and near the segment's intersections with the southern and northern caldera rim and caldera boundary faults. Source location 1 marks the surface projection of the model reservoir directly beneath the summit, whereas source location 2 marks a reservoir 5 km to the south of the summit. Data presented in Şen et al. (2003) was utilised to determine the locations of the caldera rim and major faults.