Distinct patterns of volcano deformation for hot and cold magmatic systems

Domain: earth_science
Published: 2025-01-09
Editorial outcome: major_revision

Abstract
Volcano deformation can be detected over timescales from seconds to decades, offering valuable insights for magma dynamics. However, these signals are shaped by the long-term evolution of magmatic systems, a coupling that remains poorly understood. Here we integrate thermal models of crustal-scale magmatism with thermo-mechanical simulations of ground deformation. This allows us to determine the influence of magmatic flux over 10^5–10^6 years on viscoelastic deformation spanning a 10-year observation period. Our results reveal a coupling between surface deformation and the thermal evolution of magma systems, modulated by magma flux and system lifespan. Relatively cold magma systems exhibit cycles of uplift and subsidence, while comparatively hot plumbing systems experience solely uplift. These findings align with geophysical observations from caldera systems, emphasizing the potential of surface deformation measurements as tool for deciphering the state and architecture of magmatic systems. Considering long-term magmatic system evolution is imperative for accurate interpretation of volcanic unrest.

Introduction
Analyzing ground displacement in volcanic regions is crucial for assessing volcanic hazards and enhancing our ability to forecast the initiation and progression of eruptions1,2. Space-borne Synthetic Aperture Radar Interferometry (InSAR) has revolutionized the monitoring of volcano deformation on a global scale, offering insights into diverse patterns of volcanic deformation in both space and time3–6.

However, temporal patterns of volcano deformation exhibit considerable variability, encompassing steady ground uplift over days to centuries (e.g., Laguna del Maule (Chile)7; Long Valley Caldera (USA)8; Corbetti (Ethiopia)9, and alternating uplift and subsidence cycles (e.g., Campi Flegrei in Italy10,11; Aluto (Ethiopia)12,13; Askja caldera in Iceland14; or Okmok caldera (USA)15,16). Despite rapid advances in the availability of high-quality data, understanding the underlying principles governing the vast variety of ground deformation signals at different volcanoes remains a major challenge. While phases of uplift are typically attributed to emplacement of magma intrusions3,17, a range of mechanisms to explain episodes of subsidence over timescales from months to decades have been proposed, broadly falling within four categories: (1) Pressure changes induced by fluid movement within hydrothermal systems13,18, (2) Interactions between magma pressurization and crustal heterogeneiety19, (3) Volume contraction during magma cooling, crystallization, and/or gas release20–22, and (4) Viscoelastic relaxation of either a shell surrounding a compressible magma chamber23–25 or of a viscoelastic layer below a pressurized reservoir26,27.

Viscoelastic processes have been shown to markedly impact surface deformation. Source condition estimates obtained using numerical models allowing for temperature-dependent viscoelastic rheology, require lower pressures compared to analytical solutions of elastic deformation28,29, and impact temporal trends of ground displacement on timescales of years to decades30–33. Conversely, numerical models investigating the evolution of magmatic systems over geological timescales suggest that, depending on the rate of magma supply, crustal thermal structures evolve over hundreds of thousands to millions of years34–39. Given the temperature-dependence, the viscoelastic response of rocks subjected to strain is largely governed by the long-term thermo-mechanical history of the crust40,41. While these timescales are much longer than those typically invoked in surface deformation modeling, we anticipate that magmatic systems undergo different styles of ground deformation at various stages of their lifecycle. This temporal evolution of crustal rheology may be essential in explaining the great diversity of spatial and temporal deformation behavior that is observed at different volcanoes worldwide42. Yet, even state-of-the-art thermally controlled models of volcano deformation consider the temperature field and hence viscoelastic rheology to be static, often relying on simplified geometries of magma plumbing systems25,29. This is an attempt to reconcile the vastly different timescales of ground movements and magmatic system evolution to explore these dependencies in space and time.

Here we use finite difference modeling to simulate the thermal evolution of the Earth’s crust, driven by varying long-term magmatic fluxes (see methods section for details). This is achieved through pulsed injection of magmatic sills and dikes at three different rates (high: 4.7×10^-4, mid: 3.2×10^-4, and low: 1.9×10^-4 km^3/yr magma flux cases) into the crust over a period of 1 million years (Fig. 1). Subsequently, the temperature distribution in the crust at 200,000-year intervals is integrated into a finite element code to solve for the temperature-dependent thermo-mechanical evolution, arising from a pressure source. Viscoelastic surface deformation is then monitored over a 10-year observation period. In reality, the location and size of the source would be controlled by the mechanical conditions and history-dependent stress field43,44, but for illustrative purposes, we use an oblate source situated at depths of either 5, 10, or 15 km and test a range of source sizes. Using this approach, we quantify the impact of long-term magma supply rates into the Earth’s crust over the lifespan of a crustal magmatic system (timescales ranging from 10^5 to 10^6 years) on spatial and temporal patterns of surface deformation over timescales spanning 10^0–10^1 years. Our primary focus is on the role of viscoelastic crustal rheology, reflecting dynamic changes in the evolving temperature structure of the plumbing system over time. We show that our model predictions are consistent with independent geophysical imaging of the subsurface and ground deformation time-series at two caldera systems in the East African Rift. Finally, we discuss the implications that can be drawn for volcanic deformation studies and monitoring.

Results
Magma reservoir growth
To link the long-term rheological evolution of the crust to surface deformation, we first simulated the growth and thermal maturation of magma reservoirs and their surrounding wall-rocks. Our results reveal a diverse range of thermal architectures and reservoir growth conditions primarily influenced by the long-term magma supply rate and the duration of magmatic activity (Fig. 2). Initial reservoirs accumulate at near-solidus temperatures of 700°C and gradually become hotter, reaching peak temperatures after 1 Ma of magma injection of 850–900°C (high flux case: 4.7×10^-4 km^3/yr), 800–850°C (intermediate flux case: 3.2×10^-4 km^3/yr), and 750°C (low flux case: 1.9×10^-4 km^3/yr). Generally, reservoir construction progresses from the bottom to the top over time, with initial growth occurring at a depth range approximately between 10 and 15 km, where ambient crustal temperatures are higher compared to shallower levels. In each magma flux scenario, as time progresses, the isotherms within the upper 10 km of the model domain gradually rise towards the surface, with the steepest gradient observed in the high magma flux case. Simultaneously, temperature distributions become more uniform in regions deeper than 10 km and more concentric within the growing reservoirs through time. However, strong local disturbances in the temperature field result from the last few magma pulses.

As magmatic dikes and sills inject the crust in random locations, the evolving crustal temperature field is subjected to stochastic effects. A repeated run of the intermediate magma flux case demonstrates similar reservoir growth conditions, with a persistent magma body being emplaced between 200 and 400 ka since the onset of magma injection (Figs. 2f–j; 3a–e; Fig. S1). However, the fine-scale temperature distribution is mostly influenced by the spatial history of magma injections, particularly in the earlier stages of the magmatic system.

Melt fraction plots indicate that during the early stages, the magmatic system comprises multiple pockets that coalesce into a single reservoir in later stages (Fig. 3f–j). In each case, early-stage reservoir growth is characterized by crystal mush with melt fractions <50%, increasing towards liquid-rich (i.e., eruptible) reservoir cores with about 75% melt in the high flux case, and 50% for the intermediate flux case after 1 Ma of magmatism. As shown by various modeling studies, volatile exsolution can modulate the pressurization history of magma reservoirs45,46, and deforming crystal mushes may be parameterized as poroelastic or poro-viscoelastic material47–49. Here we focus solely on thermo-viscoelastic effects considering that sensitivity analyses indicate their dominant contribution to surface deformation49 and given that the distribution and evolution of porosity in magmatic mush systems in not well understood.

The viscosity structure of the crust evolves in tandem with changes in temperature and melt fraction fields (Fig. 3k–o). Across the timeframe of magmatism from 200 ka to 1 Ma, crustal viscosities change in the low magma flux scenario from 10^25 to 10^20 Pa s at 5 km depth, 10^18 to 10^17 Pa s at 10 km, and 10^19 to 10^17 Pa s at 15 km depth (all depths at the center of the model domain). Over the same period, the intermediate magma flux case results in viscosities of 10^18 to 10^16 Pa s at 5 km depth and 10^16 to 10^15 Pa s at 10 and 15 km depth. In the high magma flux case, crustal viscosities evolve from 10^17 to 10^15 Pa s at 5 km depth, 10^16 to 10^15 Pa s at 10 km, and remain on the order of 10^15 Pa s at 15 km depth. Viscosities around 10^15 Pa s typically result in viscoelastic relaxation times (calculated as: \tau \approx \eta/(0.5*G), where \tau is the relaxation timescale in s, \eta is the viscosity in Pa s, and G is the shear modulus in Pa) of a few days or less, rendering them observable only with very high temporal resolution methods. Conversely, viscosities approximately at 10^20 Pa s correspond to relaxation timescales spanning hundreds of years, suggesting effectively elastic behavior over typical observation timescales. Viscosities falling within the intermediate range of 10^17–10^18 Pa s, however, exhibit relaxation timescales on the order of a few years, impacting our observations of volcano deformation. These considerations show that long-term magma flux and duration of magmatism exert first order controls on crustal rheologies and are therefore expected to give rise to differences in surface deformation patterns.

Long-term evolution of volcano deformation
We model short-term deformation patterns over a time span of 10 years for overpressure sources at 5, 10 or 15 km depth and at different long-term evolutionary stages of magmatic systems (i.e., from 0 to 1 Ma in 200 ka increments; cf. Figure 1). The predicted vertical displacements exhibit distinct temporal patterns influenced by the duration and magma flux of the long-term magmatism, and the depth of the overpressure source (Fig. 4). The absolute values of displacement produced in these models scale with the source overpressure applied and are overestimated by the 2-D model set-up. Similarly, cavity size has an impact on the magnitude of displacement, but not the overall pattern (Fig. S7). Thus, we focus on the temporal patterns and timescales of displacement rather than absolute values.

Ultimately, the elastic or viscous response of crustal materials to deformation is governed by the ratio of relaxation to forcing timescales (Deborah Number)43,44. Although the relaxation behavior in our model is governed by distribution of viscosities rather than single values, we do not model the injection process itself. The step-change cavity pressurization does not reflect a specific igneous process and cannot be considered equivalent to a forcing timescale due to magma injection. Instead, our study focuses on the time-dependent thermo-viscoelastic behavior of crustal rocks in response to a pre-formed, pressurized cavity.

The viscoelastic deformation history of sources at 10–15 km depth generally follows a similar trajectory: an initial rapid uplift, typically spanning 6 months to 2 years, succeeded by slower deformation. Given that magma intrusion pulses are expected to operate over timescales of days to years in most volcanic systems50, these relatively short timescales suggest that discerning surface displacement resulting from magmatic intrusion and subsequent viscoelastic deformation may pose challenges, particularly for deep sources, as also noted in previous studies29,51. This behavior remains consistent across various magma supply rates, although during the initial stages of magmatism (<400 ka), the transition between these two regimes extends over several years for low and intermediate fluxes.

In contrast, pressure sources at 5 km depth produce both uplift and subsidence (Fig. 4) with notable differences influenced by the thermal state and the rheology of the crust, which depend on the magmatic flux, the history of injection locations, and duration of magmatism (Fig. 3). The overall magnitude of surface displacement generally increases with the magma supply rate, showing peak values in scenarios featuring high magma flux, while cases with reduced magma flux result in lower displacements. For 200 ka of magma injection at low and intermediate magma supply rates, subsidence is recorded for the entire 10-year observation period, whereas, in the high magma flux scenario, initial subsidence is only observed in the first year. From 400 ka to 1 Ma, in the high magma flux scenario, uplift occurs at gradually decreasing rates. Conversely, intermediate flux scenario initially experiences a subsidence phase spanning up to 600 ka. At 1 Ma of magma injection, in the low flux scenario, subsidence is initially observed followed by uplift. Notably, the transition point from subsidence to uplift occurs progressively earlier and the total subsidence increases with longer durations of magmatism (Fig. 5).

As illustrated by a repeat run of the intermediate flux scenario (dashed lines in Figs. 4, 5b), the magnitude and temporal evolution of surface deformation can be influenced by the long-term history of injection locations, which are randomized in our modeling framework. Although the overall temporal trends are similar, the spatial injection history can significantly affect the magnitude and duration of viscoelastic subsidence and uplift patterns. More spatially focused magma injection can lead to higher local temperatures and, consequently, greater magnitudes of viscoelastic deformation (e.g., intermediate flux 800 and 1000 ka repeat runs; Fig. 5b). Therefore, the history of spatial injection locations may be an important parameter in explaining the diversity of deformation signals that occur globally.

Although, temporal patterns of surface deformation are impacted by long-term magmatic processes, coupling thermal and thermo-mechanical simulations reveals that spatial pattern of surface deformation is primarily influenced by source depth and geometry rather than the thermal architecture and growth conditions of the magmatic reservoir (Fig. S2). The initial, elastic response of the crust to the pressurized cavity is insensitive to the thermal structure, but the spatial pattern of the viscoelastic response is affected by thermal heterogeneities (Fig. 6). Thermal structures in which isotherms are irregularly distributed relative to the overpressure source results in vertical displacement patterns that develop asymmetry over time in the visco-elastic phase. In the specific cases (Fig. 6a, b), where higher temperatures prevail on the right-hand side relative to the pressure source, displacement patterns exhibit spatial differences in the extent of viscoelastic creep, evolving over time (Fig. 6c, d). This leads to asymmetry in visco-elastic displacement, relative to the elastic component of displacement at t=0 years, with a shifted center of displacement and varying magnitude response on the sides of the model domain (Fig. 6e, f). Given that the source in both scenarios remains constantly pressurized at 1 MPa, the shift can be attributed to heterogeneous viscoelastic deformation.

To summarize, viscoelastic displacement over an observation period of 10 years is influenced by the duration of magmatism, magma flux, heterogeneity of the thermal structure, and overpressure source depth. Longer durations of magmatism and higher magma fluxes generally result in greater displacement, although deeper overpressure sources exhibit similar displacement for different magma fluxes. Shallow overpressure sources show initial subsidence in some cases, particularly for shorter durations of magmatism and lower magma fluxes. This subsidence is not observed in deeper source depths. Lastly, asymmetries in the surface deformation pattern may develop over time in response to heterogeneous viscoelastic structure.

Origin of viscoelastic subsidence
The observation that several simulations with a pressure source at 5 km depth exhibit subsidence over a wide range of timescales and magnitudes is unexpected, given previous studies indicating that uplift should predominate when a constant pressure (stress) boundary load is applied28,32,33,52–54. We use the low flux case, with a pressure source at 5 km and 800 ka of magmatism to illustrate the origin of this subsidence (Fig. 7). The initial elastic deformation field at t=0 years displays a pronounced upwards component and a smaller downwards-directed vector field, resulting in net uplift. One year after the initial timestep, the downwards component has markedly increased in magnitude, caused by rapid viscoelastic creep towards the high-temperature zone underlying the shallow overpressure source (Fig. 7b). This effect leads to a net ground displacement at the surface lower than the initial elastic surface uplift, manifesting as subsidence.

Although viscoelastic creep towards the high-temperature zone beneath is still ongoing after 5 years, the process has now slowed (Fig. 7c). At this stage, the greater magnitude of the upward-directed and much slower viscoelastic creep results in ground uplift relative to the initial state.

The duration and magnitude of the subsidence period, as well as the timing of the transition from subsidence to uplift, depend on the thermal state of the plumbing system and the relative location of the overpressure source. If ambient temperatures are high, such as in the high flux case and at later magmatic stages for the intermediate flux, the effect is not observed (Fig. 4). This is because temperatures surrounding the source are high in all directions, leading to a similar viscoelastic response timescale in the upwards and downwards directions. However, if ambient temperatures surrounding the overpressure source are low and much higher in the underlying plumbing system, subsidence is observed. The duration of this process, and therefore the duration of subsidence, is thus governed by the timescale of viscoelastic creep. The occurrence of either subsidence or uplift can be predicted by analysing the difference in relaxation timescales of crustal rocks, which is primarily governed by their viscosity structure. When there is a significant difference in relaxation timescales between the rocks above and below the pressure source, subsidence is likely, as the upper layer behaves elastically while viscous relaxation occurs in the lower layer (see Fig. 8). In contrast, when the difference is minimal, both layers undergo time-dependent deformation, resulting in net uplift. Thus, the relative difference in relaxation timescales, rather than absolute values, is the key factor controlling the deformation behavior.

In this study, we do not assume any specific mechanism for pressure generation within the modeled cavity. However, two widely accepted processes for magma pressurization are the injection of fresh magma into a preexisting reservoir via a feeder dike and the accumulation of an exsolved volatile phase within the magma body50. For volatile-driven mechanisms, where pressure is generated internally, we do not expect these processes to directly influence the relaxation behavior of the crustal rocks. However, magma transport through a feeder dike system could affect the stress field beneath the intrusion, potentially counteracting the downward-directed flow of hot host rocks. More complex models are needed to explore the interaction between magma injection and viscoelastic crustal rheology to determine whether this hypothetical effect would be significant. Nonetheless, petrologically constrained magma recharge timescales, typically ranging from days to months50, are shorter than the 10-year observation window modeled in this study. Therefore, our results are applicable to scenarios in which the injection process has terminated.

Our results are broadly consistent with the mechanism suggested by Yamasaki et al.26 to explain post-elastic subsidence, namely, downwards-directed viscoelastic creep resulting from a heterogeneous thermo-mechanical architecture of the crust. However, while Yamasaki et al.26 focus on a two-layer model where injections are confined to the boundary between an elastic and viscoelastic crustal layer, we account for different rheological architectures of the crust resulting from different magma supply rates and durations of magmatic activity. This allows to map out under which conditions of reservoir growth and magma system lifespan such effects are expected to prevail. Downwards-directed viscoelastic creep is distinct from subsidence mechanisms found in viscoelastic shell models23–25, which show subsidence only in calculations with high magma compressibility and short magma injection durations compared to the viscoelastic relaxation timescale. These models do not account for viscoelastic creep in different spatial dimensions. Our results show that post-elastic subsidence can occur independent of the duration of the injection episode and independent of magma compressibility. However, in accordance with the mechanism proposed by Townsend25, post-elastic subsidence does not necessarily involve a volume change of the magma reservoir.

Comparison to geophysical observations
While the particular mechanism generating surface deformation is often difficult to discern at volcanoes and not seldom fiercely debated, our model produces testable predictions on the relationship of temporal surface deformation patterns and the thermal architecture of magmatic feeding systems that can be verified against real-world observations. Our results define distinct visco-elastic phases of subsidence followed by uplift in relatively cold magmatic systems, contrasting with hot systems exhibiting a singular pattern of ground uplift (see Figs. 4, 5; Table S1). We test these model predictions against the findings of recent geophysical imaging of the igneous feeding systems and the observed surface deformation patterns at two neighboring East African volcanoes: The Aluto and Tulu Moye volcanic complexes.

Both Aluto and Tulu Moye are large, nested caldera complexes, located approximately 50 km apart in the Central Main East African Rift55–57. They share several prominent features, such as bimodal basalt-rhyolite magma suites, giving rise to large magnitude explosive eruptions and smaller effusive events13,58. Both volcanoes have extensive eruptive histories, spanning at >300 ka for Aluto13 and >100 ka for Tulu Moye59. The most recent activity at Aluto, marked by an explosive event, occurred about 400 years ago13 and historical accounts and eyewitness reports suggest emplacement of a rhyolitic lava flow at Tulu Moye around the year 1900 CE57,60. Despite these similarities, recent 3D magnetotelluric imaging has revealed contrasting states of the magmatic plumbing system beneath the two volcanoes61.

The magnetotelluric method images the electrical conductivity structure of the subsurface based on natural electromagnetic field variations. Anomalies of higher electrical conductivity beneath volcanoes can be attributed to variations in temperature, the presence of melt, saline geothermal fluids, and/or hydrothermal alteration minerals62,63. At Aluto and Tulu Moye, an overall similar conductivity structure has been imaged (Fig. 9; refs. 61,64,65.). This trifold architectural pattern, consisting of a shallow-level anomaly interpreted as a hydrothermal alteration cap, underlaid by intermediate and deep-level conductivity anomalies, is a distinctive feature observed in numerous volcanoes worldwide63. However, at Aluto and Tulu Moye, significant differences exist in the magnitude of the intermediate (~4–6 km depth) and deep (>8 km depth) conductivity anomalies, showing higher conductivity values for Tulu Moye.

Combining petrological phase equilibria and electrical conductivity modeling, these differences in conductivity can be attributed to temperature and melt fraction differences for the anomalies beneath each volcano61. At Aluto, temperatures of 725–745°C have been inferred, while modeling the magnetotelluric signal at Tulu Moye indicate temperatures between 785–905°C. As shown in Fig. 10 and Table S2, these temperature ranges are well represented by the temperature distributions of our simulations. We suggest that average temperatures of the magma above solidus are most relevant to the comparison given transient effects that impact extreme values, as well as considering the resolution of magnetotelluric imaging on spatial scales of ~10^2–10^3 m and the consequent averaging of the temperature distributions. Generally, model runs with average temperatures within the range observed at Aluto show either subsidence or subsidence that transitions into uplift at a later stage, while simulations with mean temperatures corresponding to Tulu Moye show monotonous uplift (Fig. 10; Table S2; Fig. 5). Matching of our predictions and temperature distributions for Aluto and Tulu Moye thus suggests that these differences could be driven by either contrasting durations of magmatic activity for the two volcanoes at the same long-term magma input rate (e.g., high flux scenario 200 ka versus 800 ka; Fig. 10; Table S2), or reflect magmatic systems at a similar evolutionary stage that have been constructed with contrasting long-term magma addition rates (e.g., low flux scenario versus high flux case at 600 ka). Long-term crustal magma flux estimates at Tulu Moye and Aluto are, however, not available and published geochronological data does not permit us to distinguish between the different possibilities.

Surface deformation has been detected using InSAR at both Tulu Moye and Aluto (Biggs et al. 2011). Between 2004 and 2010, Aluto showed two cycles of uplift lasting several months, each followed by longer periods of subsidence (Fig. 9). In contrast, no subsidence has been observed at Tulu Moye, but solely ground uplift has in the period between 2008 and 2010, as well as since 2016 (refs. 12,66; COMET Volcano Deformation Portal: https://comet.nerc.ac.uk/comet-volcano-portal/; Fig. 9). Estimates of the deformation source depth and geometry for both volcanoes were carried out using elastic analytical solutions. The deformation pattern at Tulu Moye can be explained by a sill-like geometry at about 7.7 km depth beneath the surface67 and at Aluto by a point source at a depth range between 4.8 and 5.4 km13. These observations are fully consistent with our modeling results suggesting continuous uplift pattern for hot plumbing systems like Tulu Moye and subsidence-uplift cycles for colder magmatic systems like Aluto (Figs. 8, 9; Table S2). Comparing the deformation source estimates to the thermal structure of the plumbing system beneath the two volcanoes reveal that both pressure sources are underlain by a high temperature zone but also that the ambient temperatures at the level of the pressure source are higher in the case of Tulu Moye. In Fig. 9, we present two simulations: one for Aluto (intermediate flux, 200 ka) and another for Tulu Moye (high flux, 800 ka), both matching the observed temperature distributions (Fig. 10). Although multiple simulations could fit the temperature profiles for each volcano, the cases matching Tulu Moye consistently show monotonous uplift, leading to similar deformation time series. For Aluto, we selected a simulation depicting subsidence over timescales comparable to the observed deformation series, acknowledging that other simulations might show a higher frequency of uplift-subsidence cycles compared to what we observe at Aluto. Although our results do not rule out any contribution of hydrothermal processes13 in shaping the deformation, we note that the source depth >5–8 km are difficult to reconcile with the imaged extent of hydrothermal alteration <2 km. We therefore propose that the differences in deformation regime for Tullo Moye and Aluto are primarily a reflection of the contrasting thermo-mechanical architecture of the magmatic plumbing systems.

Implications for volcanology
The findings presented in this study carry significant practical implications for understanding volcanic unrest signals. Our simulations indicate that spatial deformation patterns can undergo time-dependent shifts when magma bodies interact with a heterogeneous thermal structure of the crust. Misinterpreting such changes in spatial deformation patterns resulting from viscoelastic creep of crustal rocks as indications of underground magma movement could result in erroneous forecasts. Furthermore, our model calculations reveal characteristic temporal patterns of surface deformation for cold and hot magmatic plumbing systems, validated by fully independent geophysical imaging and InSAR data. This implies that deformation cycles, which are widely observed at caldera volcanoes, may mirror the underlying thermo-mechanical architecture of the plumbing system. This is important as it offers the opportunity to use differences in surface deformation time series to gain insights into large-scale magmatic system architecture, going beyond conventional deformation source inversion.

Recent studies highlight the vital role of viscoelastic processes in volcano deformation25,27,29,32,43,44,68, and our findings stress the importance of considering these processes more broadly, especially in shallow and young magmatic systems. The thermal architecture of transcrustal magmatic systems is formed by the repeated injection of numerous magma batches over prolonged periods of time. Unlike viscoelastic shell models or geometries with stacked elastic and viscoelastic layers, such systems exhibit a continuum of elastic and viscoelastic materials. This suggests that even shallow deformation sources (e.g., <5 km depth), often assumed to behave elastically, may be significantly impacted by time-dependent rheology with heterogeneous expression in different spatial dimensions. Incorporating temperature-dependent viscoelastic processes into the geodetic inversion of volcanic deformation signals is crucial for accurate monitoring. Neglecting these factors could lead to inaccuracies in source depth or volume calculations, highlighting the necessity of integrating them into analysis methods. Such integration could effectively be based on the results of geophysical imaging, such as magnetotelluric and seismic methods, providing information on the distribution of crustal temperatures beneath volcanic systems.

In addition to these practical aspects, our study holds further implications for fundamental volcanology. We show that the architecture of magmatic plumbing systems, as governed by the long-term magmatic flux, the history of injection locations, and system lifespan. Given our model prediction that monotonous uplift and subsidence-uplift cycles correspond to hot and cold magmatic systems respectively, these differences could be used for first-order estimates of subsurface magmatic conditions, such as temperature distributions and potentially long-term magma fluxes, which are otherwise difficult to estimate.

A major challenge in this conceptual framework is that multiple trade-offs exist between different variables that impact surface deformation fields. However, estimates of the lifespan of magmatic systems, magma volatile contents, depth range of injection, and long-term magma input rates can be recovered through geochronological and petrological methods69–71. Integrating age determinations and petrological studies with thermal and deformation modeling has therefore great potential to improve our understanding of the architecture and dynamics of igneous plumbing systems. This is crucial because, while surface deformation is now routinely measured on a global scale, our understanding of crustal magma budgets and the thermal state of volcanic plumbing systems remains limited.

In summary, our findings underscore how the evolving thermal and rheological architecture of magmatic systems over their lifecycle shape short-term volcano deformation dynamics, with far-reaching implications for fundamental volcanology, hazard mitigation and volcano monitoring strategies.

Methods
Modeling long-term magma evolution
Numerical simulations were conducted to investigate the long-term thermal evolution of magmatic systems and the surrounding wall rocks. The simulations involved the pulsed injection of dikes and sills and were ran using the open-source software package ‘MagmaThermoKinematics’ (MTK), implemented in the high-level programming language Julia (Schmitt et al. 2023). A comprehensive description of the software, including its source code and benchmarking against previously published thermal models69 of pulsed magma injection, is provided in ref. 70 and is available at https://github.com/boriskaus/MagmaThermoKinematics.jl. Here we summarize the general methodology, describe the setup and parameters used in our specific case.

The core of the simulation is based on solving the heat diffusion-advection equation using the finite difference method. The governing equation used in Magma-thermokinematics can be expressed as:

\rho c_p \frac{\partial T}{\partial t} + v_j \frac{\partial T}{\partial x_j} = \frac{\partial}{\partial x_i} k \frac{\partial T}{\partial x_i} + H_r + \rho Q_L \frac{\partial \theta_s}{\partial t} (1)

Here, T represents temperature in Kelvin, \upsilon denotes the host-rock velocity in response to the injection of new sills or dikes, t signifies time in seconds, \theta_s stands for the solid fraction, \rho is the density, c_p denotes heat capacity, k represents the thermal conductivity, H_r is the radioactive heat source, and Q_L is the latent heat. When magma intrudes, it displaces the host rock to make room for itself. To model this process, an injection routine based on an analytical solution is used that inserts the dike temperature into the temperature field and defines the velocity field needed to create a penny-shaped crack in an elastic half space. To address numerical stability issues arising from the time derivative of the solid fraction (\theta_s), MTK locally adjusts the heat capacity to incorporate latent heat release, thereby enhancing the stability of the solution. Furthermore, since material properties, such as melt fraction, are temperature-dependent, a non-linear iteration loop is employed.

The software discretizes the time derivative explicitly, using the following stable timestep criterion:

\Delta t < \frac{1}{10} \frac{\Delta x^2, \Delta y^2}{\kappa} (2)

where \Delta x^2 and \Delta y^2 are the grid spacings, and \kappa is the thermal diffusivity (\kappa=k/(\rho c_p)). A staggered grid finite difference scheme has been employed to compute the spatial derivative, effectively addressing variations in thermal conductivity between temperature grid points. The implementation incorporates a temperature-dependent thermal conductivity72.

Simulations were carried out in 2D geometry to track the spatial development of thermal heterogeneity within the crust. Two sets of geometries were employed to assess the impact of different domain sizes: (1) A geometry with a physical domain size of 20×20 km, discretized with 300×300 numerical nodes. (2) A larger geometry with an expanded horizontal physical domain size of 40 km and a depth of 20 km, discretized by 600×300 numerical nodes. Initial testing showed that these geometries provide similar results (Fig. S4). We thus focus in the following only on the smaller geometry case.

A generic melting parameterization, implemented in MTK based on ref. 71 was used for this study. In the used parameterization, melt fraction is calculated from the temperature in °C as:

\psi = \frac{(800-T)}{30}

\theta_{melt} = \frac{1}{(1+e^\psi)}

with a liquidus-solidus range of 1000–700°C, which is applicable for intermediate to silicic magma compositions. The initial temperature profile was determined by applying a linear geothermal gradient of 25°C per kilometer. Flux-free boundary conditions were applied to the sides of the model. To explore various magma injection rate, we adjusted the fluxes by modifying the periodicity of sill and dike emplacement. In this study, we focus on three distinct scenarios labeled as follows: high magma flux case (4.7×10^-4 km^3/yr), intermediate flux (3.2×10^-4 km^3/yr), and low magma flux case (1.9×10^-4 km^3/yr). These magma fluxes have been chosen to capture a wide range of crustal thermal histories with incubation periods for magma reservoir growth ranging between ~200 ka (high flux case) to ~1000 ka (low flux case) (Fig. 2). The volume of dikes and sills was computed based on the width (200 m) and length (1500 m) of injected dikes, assuming a penny-shaped geometry. Consequently, individual injections resulted in a volume of 0.94 km^3. Injection locations were randomly distributed over a height of 10 km and a width of 6 km relative to the center of the entire domain. Injection angles of dikes and sills were randomly varied from 80 to 100 and −10 to 10°, respectively. To assess the impact of random injection on the thermal evolution of the magmatic system, the intermediate flux scenario was repeated using a different randomized set of injection locations. All simulations ran over a total duration of 1 million years. The evolving temperature field was tracked at intervals of 200,000 years and subsequently utilized as input for crustal deformation modeling.

Surface deformation modeling
To compute short-term surface deformation pattern at different stages of the evolutionary history of magmatic systems, we coupled the thermal output with thermo-mechanical models. Crustal deformation modeling was done using the structural mechanics and heat transfer modules of the finite element code COMSOL Multiphysics version 5.1. The model setup is based on the benchmark procedure for surface deformation developed in ref. 53. We used a 2D geometry, in order to determine the impact of spatial heterogeneity of the temperature field on the results, which is in not easily accounted for in the more commonly used axisymmetric model formulations. We note that the 2D geometry overestimates the magnitude of deformation compared to previous benchmark models53, but this is less problematic as we are interested in general trends and differences between simulations at various evolutionary stages of the magmatic system rather than modeling deformation at a specific volcanic system.

Model geometries of either a 20×20 km or 20×40 km were used equivalent to the thermal simulations. First, the temperature field corresponding to a particular evolutionary stage from the thermal model was incorporated using an interpolation function. The geometry was then parameterized as a temperature-dependent viscoelastic rheology using a Standard Linear Solid material model with the following properties: Heat capacity of 1050 J/kg K, density of 2700 kg m^-3, thermal conductivity of 2 W/(m K), Bulk modulus of 13.3×10^9 N/m^2, and shear modulus of 8×10^9 N/m^2. In addition, a subset of calculations was carried out with generalized Maxwell rheology, using the same set of parameters (Fig. S5). For reference, a comprehensive list of parameters employed in the thermal and thermo-mechanical simulations is given in Table S1. Additionally, the modeling setup is illustrated in Fig. S5.

Viscosity was calculated using an Arrhenius law of the form:

\eta = A_d e^{(E_a / (R A T_e))}

where A_d is the Dorn parameter of 10^9 Pa s, E_a is the activation energy of 120,000 J/mol, R is the gas constant of 8.134 J/(mol K), and T is the temperature field from the thermal simulations. The relaxation time was calculated as:

\tau = \frac{\eta}{G_{inst} \mu}

where \eta is the viscosity, G_{inst} is the instantaneous shear modulus, and \mu is the fraction (0.5) of the shear modulus of the viscoelastic arm.

Overpressure sources were modeled as ellipsoidal cavities with dimensions of 800×200 m at a depth of either 5, 10 or 15 km, unto which a constant boundary load of 1 MPa was applied if not stated otherwise (e.g., Figure 9; Figure. S6; Fig. S3). Given the impact of cavity size on surface deformation (e.g., refs. 43,44), we also explored smaller and larger source dimensions as shown in Fig. S7. Conceptually, we envision that the modeled overpressure sources represent either single intrusive bodies (e.g., shallow laccoliths), pressurized melt lenses in a crystal mush, or recharge events into a mush zone without invoking magma mixing or mingling processes. It is important to note, however, that we do not explicitly model the intrusion process that lead to cavity formation. Consequently, the initial elastic response should not be interpreted as resulting from a specific igneous process. Instead, our focus is solely on the thermo-viscoelastic processes resulting from a pre-formed, overpressurized cavity. The surface was modeled as a free boundary, while roller boundaries were applied to the sides of the model domain, and a fixed constraint corresponding to the initial condition was applied to the bottom of the model domain. The model geometry was discretized with typically >31,000 triangular mesh elements.

We solved the temperature-dependent crustal deformation field in two steps following ref. 53. In a first step, the temperature field from the thermal modeling results, corresponding to a particular evolutionary stage (e.g., 400 ka of magmatism), was implemented through the heat transfer module in COMSOL Multiphysics. In a second step, crustal deformation in response to the overpressure source was solved for a duration of 10 years with step size of 6 month using the solid mechanics solver. We assume that the temperature field can be treated as quasi-stationary over the timespan of interest (10 years). The COMSOL.mph and thermal input files are available in a zenodo repository.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In their paper \"Distinct patterns of volcano deformation for hot and cold magmatic systems\" Wenner et al. develop a long term (1 Ma) thermal model of crustal magma system establishment (using low, mid and high background magmatic fluxes) and evolution, based on random dike and sill injections into a predefined 2D region to assess temperature, melt-fraction and viscosity of the medium over time. These results are used as background properties into which ellipsoidal pressure sources and 5,10,15 km are embedded (different model, COMSOL), to model (vertical) displacements due to an impulse pressurization and the following visco-elastic response of this heterogeneous medium observed over a 10-year timespan. Depending on stage in the evolution (200-1000 ka) and the resulting system temperature that also depends on magmatic flux, the authors find substantially different deformation responses that range from expected post-eruptive exponentially decaying uplift for hot systems to subsidence and uplift episodes for colder systems. To explain this, they invoke viscous creep and its interaction with (a) an evenly hot surrounding material, which gives an equal viscoelastic response in all directions, and (b) a system that is hotter below than above, resulting in downwards viscous creep first and subsidence. They are able to reproduce the deformation behaviour at two systems in the East African Rift (Aluto and Tullu Moye) using MT observations to constrain the background temperatures and inferences about the evolution stages of the systems / long term magma flux.\nThis paper puts forward an interesting mechanism to furhter our understanding of magmatic systems and the sometimes puzzling deformation fields we are concerned with. They rightfully point out that they report a behavior that has not been captured in analytical models and support through their results the importance of considering heterogeneously warmed crust and the importance of where injections occur inside the background temperature field. I guess a main issue with their model is the dependence on both background temperature and long-term magmatic flux, which are both difficult to constrain and result in overlapping scenarios that can create different deformation fields (see their East Africa example).\nI believe, this concept is important to convey to the broader community and will be of great interest, which warrants publication in Nature Communications in my view. My comments below probably place the necessary revisions somewhere between minor and major, somewhat depending on what the figures show etc.\n1) There's a lot to consider in this paper: magma flux, evolotionary stage of the magmatic system, depth of the injection, or generally location of injection into the background temperature field. Clearly, a short form paper and touch on only so many of these aspects, but I believe the paper would benefit from some reconsiderations of cases that are considered: It seems that endmember cases might be best suited to convey the concepts: high and low flux, shallow and deep injections (and perhaps furhter consideration of inside the hot area and adjacent to it). As the paper is written, it is hard to figure out which change results in what deformation field. Figure 6 is a good example: the authors change the location of the injection, but also the background temperature field AND the evolution stage of the system AND the background flux. They title this figure \"Influence of thermal heterogeneity on spatial surface displacement\" ... it seems this might be well achieved by using the same, somewhat asymmetric background temperature field and injecting centrally and adjacent to demonstrate its importance relative to the location of the magma recharge. It might well be that I am missing something here, that in itself might be cause to restructure a bit.\n2) To this point - Figure 4 shows the importance of the injection history (dotted line) on the deformation response. I assume this trades off with the injection location? While I agree that the dotted line shows similar behavior for the intermediate case, it seems that the subsidence nature is lost sooner (compare solid and dotted orange lines for 5 km case for 800 and 1000 ka). This brings the same flux closer to the high flux scenario, solely based on injection location, I presume, further complicating the tradeoffs between parameters mentioned above. It seems the authors should point this out more explicitely than is currently done in L135-140, which doesn't mention the later stages of magma system evolution.\n3) If the authors retain the intermediate flux case, it should be retained in Figure 5, too.\n4) Figure 7: it would be good to show the two contrasting cases here (high/low flux).\n5) For the East Africa example, I am not able to follow why they chose the respective scenarios. The text is somewhat vague and does not specify clearly (maybe I couldn't find it) why the 200 and 600 ka and cold / hot scenarios are chosen (shown in figure 8) instead of the low and high flux after 1 Ma, or intermediate flux after 400 ka and 1Ma, which are mentioned in the text. Perhaps best would be a comparison of deformation predictions based on several models to demonstrate the impact of these choices (but maybe the deformation histories will not be matched as well?). Presumably, Aluro went through 2 injections during that time period (Y1, Y4)? Can't find that specified in the text. I would also be really curious if they've found a data set that shows the migration of the uplift center as suggested in figure 6b,d,f.\n6) I am bit puzzled by the limitations on \"implications for volcano monitoring.\" While this is an important aspect, it seems that this work also has important implications on basic volcano science. Could the authors show how we could use different deformation fields to infer heat / melt fraction distributions (even to first order), or long-term magmatic flux for these systems? That is, can we learn new things beyond the classic elastic homogeneous isotropic half-space knowledge we commonly get from deformation fields with their insights? If so (and I think we can, but there are substantial trade-offs, see above), this should be pointed out in some conceptual sense as it would increase the paper's impact.\nMinor issues:\nLine 9-10: why just forecasting?\nIntro: might perhaps be worth including non InSAR references, given the long timeseries that we are accumulating at some GNSS monitored volcanoes, particulalry highlighting these processes?\nL48: numerical model_s_\nequation 2: is there a < missing?\n503: was -> were\nsometimes Pas or Pa s ... use the latter\nL158-161: could be observable with GNSS (or tilt) if we knew what we're looking for...\nFigure 4: dotted orange line really hard to see...\nL209: add dot after \"rates\"\nL233 - 238: should be observable in data - seen anywhere?", "reviewer_2": "This manuscript presents interesting results that highlight the role of thermos-elastic deformation. I find of particular noteworthiness the results around the influence of thermal heterogeneity and how this produced corresponding heterogeneous viscoelastic deformation responses. The work will be of significance to anyone in the field of volcano deformation.\nI see no major flaws in the method or analysis, and comments/suggestions for improvement can be found in the attached file. These are minor and involve three main themes.\n1. A “toning” down of some of the claims in the introduction. This work included novel aspects, but there have been several studies with comparable elements that have looked at non-static thermally controlled deformation models of volcanic systems (a good few of which ae cited by the authors).\n2. I got there in the end, but I did find figure 4 initially difficult to understand, I may have not read thoroughly enough the opening few pages, or not realised the importance, but I think this needs work to make sure thar the results are supported by the description of the methodology, a really quick win would be referring back to Figure 1 when discussing Figure 4.\n3. I have no issue with the underlying methodology, it is appropriate and understandable. I would however like to see more commentary on the creation of the 800m x 200m pressurised cavity, or specifically more commentary on what it is not. I may be being oversensitive to this, but the manuscript is presented as b the results being envisage as created from a distinct intrusion or equivalent, after a given period of magmatic activity. However, the result does not really model the intrusion of anything, as the cavity is pre-formed in the model mesh and then pressurised. Any comments therefore on the initial elastic response cannot be overly interpreted, and the appropriateness of if representing an intrusion is also limited. The strength, as mentioned previously is looking at the viscoelastic response, and the authors already allude to this in lines 147-148. I would like to see it more explicitly expressed that the initial elastic response does not represent any expected volcanic process, while emphasising the observations on the thermoelastic response.\nThese are minor points and mainly clarifications, I have made a few suggestions around cleaning up the text in the introduction, otherwise I enjoyed reading the manuscript which details an interesting and useful piece of work.", "reviewer_3": "This manuscript takes a creative approach to the problem of reconciling different timescales involved in volcano deformation. The authors model 100s kyr thermal evolution in a crustal section, which then provides initial conditions and defines static but spatially heterogeneous material properties in a 1-10 year viscoelastic model of a single spheroidal intrusion. The approach merges notions of thermal conditioning and long-term intrusive flux with the cavity-type descriptions of magma chambers used to model active volcano deformation.\nThe numerical modeling is 2D, a mash-up of a commercial, non-open source FEM code (COMSOL) and an open source finite difference code. The authors identify some interesting and somewhat non-intuitive results that suggest a role for thermal maturity - and by extension the history of magmatic flux - in volcano deformation timeseries, including scenarios for which intrusion results in subsidence due to viscoelastic effects and cyclic behavior. The authors apply these results in a qualitative way to data from volcanoes in the East African Rift.\nI support publication – depending how the authors address the comment listed in the first below. But overall, it would be nice to see a more thorough parameter exploration to argue the main points convincingly. The authors could also do a better job justifying various assumptions in the model, as they are not unique and future work will benefit from a clearer exposition.\nBefore I believe the results, I’d like to know whether the subsidence you observe in models is simply a result of inconsistent model initial conditions.\nYou are modeling long-term assembly in a kinematic way, but the system has mechanical memory and stresses associated with emplacement of intrusions have long-range influence. Rucker et al., (2022) found that steady state cavity size varied with forcing period because viscoelastic creep is frequency dependent. This implies that the “size” of a magma chamber can’t simply be imposed kinematically in time-dependent viscoelastic problems… it’s a model result. My impression is that most viscoelastic chamber models do not address mechanical initial conditions. But chamber volume impacts surface deformation, so this impacts the significance of your results. Similar issues relating to pre-stress have been long recognized in earthquake models (e.g., Day, 1982). We’re sortof behind the curve in volcanology.\nI’m also curious whether you’d see subsidence even in the absence of sill emplacement… since you’re coupling two different codes (two entirely different numerical methods!), self-consistency must be demonstrated. It could be that one needs to derive initial equilibrium mechanical conditions in an iterative manner (similar to what is done with gravitational loading from topography in elastic problems, e.g., Cianetti et al., 2012).\nOther comments:\nYou use Standard Linear Solid rheology but it’s a little unclear how you’re deriving a single relaxation timescale (e.g., line 159). I am actually a bit skeptical of the need to use SLS rheology here despite several recent papers that advocate for it. SLS seems more a way to get sensible deformations from a physically implausible model setup (constant step chamber pressure increase implies constant magma flux into the chamber, does it not? Not even an oversimplified magmastatic approach to mass influx implies that…) rather than anything about crustal rheology. Spatially variable Arrhenius temperature dependence and time variable forcing in a Maxwell model will trade off with band-limited effects from the SLS model. Shouldn’t we use the simplest model until real data dictate otherwise? Regardless, this all needs to be tightened up, explicit rheologic and mass balance assumptions stated etc.\nAnd its not simply the relaxation time that governs deformation in viscoelastic systems as described on line 159. It’s the ratio of relaxation time to forcing timescale (the Deborah number, for Maxwell rheology) that determines whether the surface deformation reflects elastic or viscous response (Rucker et al., 2022; Liao et al. 2023). I think you could do a more convincing job with this problem if you developed the appropriate nondimensional numbers and presented your results in this way. For example, it would be great to say something more concrete about how the spatial structure of the thermal field impacts the mechanical evolution (paragraph starting ln 223). What sets the lateral extent of the chamber or the thermal preconditioning here? It plays a big role. Then I’d seek to explicitly tie the influx rate to the thermal preconditioning and then to the deformation. The latter is our main observation, and it would be nice to know whether anything unique is encoded there. Could be the answer is no.\nIts worth noting that Liao et al (2023) has essentially resolved the dichotomy between older “shell” models and spatially variable thermoviscoelastic models of magma chamber deformation. The “shell” can be precisely defined and derived analytically for each frequency component of forcing, and one can associate the “shell” spatial extent with a Deborah number contour (even for step function pressure as modeled here). That’s a fairly big result in my opinion.\nFinally, perhaps this is more of a philosophical comment, but I’m struck by some phenomenologic similarities between this manuscript and that in Karlstrom et al., (2010). For example, Fig 4 in that paper doesn’t look so different to me than Fig 3 in the present manuscript (appropriately blurring the eyes). Certainly the 2D stochastic approach to long-term recharge and notions of thermal pre-conditioning as governing short-term dynamics are similar. There are so many benefits to the modeling framework developed here of course, but I suspect that the basic phenomenology of your results could be (and to some extent already have been) more simply demonstrated.\nGreat, provocative work regardless.\nLeif Karlstrom"}, "author_response": "We appreciate the constructive suggestions regarding the clarity of our manuscript and\nmade several revisions to address these points.\nFirstly, we have revised Fig. 6 to better illustrate the impact of thermal heterogeneity\nwhile keeping other variables constant. In the revised version, we show two scenarios\nwhere injections are made centrally and adjacent to the hot area, using the same long-\nterm magma flux and evolutionary stage of the system.\nSecondly, we have revised the text to improve clarity regarding the impacts of varying\ninjection locations, background temperature fields, evolutionary stages of the system,\nand magma fluxes. Specifically, we have updated the discussion in Lines 269-277 to\nprovide a more explicit explanation of how each variable influences the resulting\ndeformation field.\nGiven the feedback from other reviewers, particularly Reviewer 3's request for an\nexpanded parameter space exploration, we have retained the intermediate flux scenario.\nHowever, we have moved any additional parameter explorations beyond those directly\nrelevant to the central findings into the supplementary materials.\n2) To this point - Figure 4 shows the importance of the injection history (dotted line) on the\ndeformation response. I assume this trades off with the injection location? While I agree that\nthe dotted line shows similar behavior for the intermediate case, it seems that the subsidence\nnature is lost sooner (compare solid and dotted orange lines for 5 km case for 800 and 1000\nka). This brings the same flux closer to the high flux scenario, solely based on injection\nlocation, I presume, further complicating the tradeoffs between parameters mentioned above.\nIt seems the authors should point this out more explicitely than is currently done in L135-140,\nwhich doesn't mention the later stages of magma system evolution.\nWe agree that the effect of injection site history on surface deformation could be more\nclearly highlighted in the manuscript. The relevant section of the text has been revised\nand expanded:\n“As illustrated by a repeat run of the intermediate flux scenario (dashed lines in Figs. 4,\n5b), the magnitude and temporal evolution of surface deformation can be influenced by\nthe long-term history of injection locations, which are randomized in our modelling\nframework. Although the overall temporal trends are similar, the spatial injection\nhistory can significantly affect the magnitude and duration of viscoelastic subsidence\nand uplift patterns. More spatially focused magma injection can lead to higher local\ntemperatures and, consequently, greater magnitudes of viscoelastic deformation (e.g.,\nintermediate flux 800 and 1000 ka repeat runs; Fig. 5b). Therefore, the history of spatial\ninjection locations may be an important parameter in explaining the diversity of\ndeformation signals that occur globally.” Lines 234-242.\n3) If the authors retain the intermediate flux case, it should be retained in Figure 5, too.\nThe intermediate flux scenario has been added to the figure.\n4) Figure 7: it would be good to show the two contrasting cases here (high/low flux).\nThe figure has been modified as suggested.\n5) For the East Africa example, I am not able to follow why they chose the respective\nscenarios. The text is somewhat vague and does not specify clearly (maybe I couldn't find it)\nwhy the 200 and 600 ka and cold / hot scenarios are chosen (shown in figure 8) instead of the\nlow and high flux after 1 Ma, or intermediate flux after 400 ka and 1Ma, which are\nmentioned in the text. Perhaps best would be a comparison of deformation predictions based on\nseveral models to demonstrate the impact of these choices (but maybe the deformation\nhistories will not be matched as well?). Presumably, Aluro went through 2 injections during\nthat time period (Y1, Y4)? Can't find that specified in the text. I would also be really curious\nif they've found a data set that shows the migration of the uplift center as suggested in figure\n6b,d,f.\nThe manuscript was not sufficiently clear in pointing out why we chose specific\nsimulations for comparison with the natural data. We now state in detail that this\ncomparison is based on matching of temperature distributions in the model and\npublished temperature estimates for the two volcanoes (Lines 382-397). To ease the\ncomparison, a new Fig. 9 and Table S1 has been added showing temperature\ndistributions and summary statistics of modelled magmatic temperature distributions in\ncomparison to the published estimates for Aluto and Tullu Moye.\nThe specific choice of simulations which are shown in Fig. 8 is now discussed is Lines\n441-448:\n“In Fig. 8, we present two simulations: one for Aluto (intermediate flux, 200 ka) and\nanother for Tullu Moye (high flux, 800 ka), both matching the observed temperature\ndistributions (Fig. 9). Although multiple simulations could fit the temperature profiles\nfor each volcano, the cases matching Tullu Moye consistently show monotonous uplift,\nleading to similar deformation time series. For Aluto, we selected a simulation depicting\nsubsidence over timescales comparable to the observed deformation series,\nacknowledging that other simulations might show a higher frequency of uplift-\nsubsidence cycles compared to what we observe at Aluto.”\nWith regards to modelled injections for Aluto: We mention in Line 411 that Aluto is\nmodelled with 2 consecutive injection pulses. These are also shown in supplementary\nFig. 3.\nRegarding the migration of the uplift centre, we have checked at some possible\ncandidates and see hints of a migration in the dataset at Tullu Moye. However, the\nquality of the data is not good enough to be included here - there are data gaps due to\nvegetation, which means the results are not conclusive. We will keep looking and if we\ndo find something, it would be an excellent topic for a follow-up paper.\n6) I am bit puzzled by the limitations on \"implications for volcano monitoring.\" While this is\nan important aspect, it seems that this work also has important implications on basic volcano\nscience. Could the authors show how we could use different deformation fields to infer heat /\nmelt fraction distributions (even to first order), or long-term magmatic flux for these systems?\nThat is, can we learn new things beyond the classic elastic homogeneous isotropic half-space\nknowledge we commonly get from deformation fields with their insights? If so (and I think we\ncan, but there are substantial trade-offs, see above), this should be pointed out in some\nconceptual sense as it would increase the paper's impact.\nThank you for this remark. A conceptual discussion of how the findings of our study\ncould be used to infer magmatic conditions has been added:\n“In addition to these practical aspects, our study holds further implications for\nfundamental volcano science. We show that volcanic surface deformation time series\ncan reflect thermal and rheological architecture of magmatic plumbing systems, as\ngoverned by the long-term magmatic flux, the history of injection locations, and system\nlifespan. Given our model prediction that monotonous uplift and subsidence-uplift\ncycles correspond to hot and cold magmatic systems respectively, these differences\ncould be used for first-order estimates of subsurface magmatic conditions, such as\ntemperature distributions and potentially long-term magma fluxes, which are otherwise\ndifficult to estimate.\nA major challenge in this conceptual framework is that multiple trade-offs exist\nbetween different variables that impact surface deformation fields. However, estimates\nof the lifespan of magmatic systems, magma volatile contents, depth range of injection,\nand long-term magma input rates can be recovered through geochronological and\npetrological methods64,65,67. Integrating age determinations and petrological studies\nwith thermal and deformation modeling has therefore great potential to improve our\nunderstanding of the architecture and dynamics of igneous plumbing systems. This is\ncrucial because, while surface deformation is now routinely measured on a global scale,\nour understanding of crustal magma budgets and the thermal state of volcanic\nplumbing systems remains limited.” Lines 484-500\nMinor issues:\nLine 9-10: why just forecasting?\nAdded “… magma dynamics …”.\nIntro: might perhaps be worth including non InSAR references, given the long timeseries that\nwe are accumulating at some GNSS monitored volcanoes, particulalry highlighting these\nprocesses?\nWe did include a wide range of GNSS studies in the introduction, including references\n8, 10-11, 14-15.\nL48: numerical model_s_\nThanks for spotting this. Corrected.\nequation 2: is there a < missing?\nYes, this was missing. It has been added.\n503: was -> were\nCorrected.\nsometimes Pas or Pa s ... use the latter\nThis has been corrected throughout the manuscript.\nL158-161: could be observable with GNSS (or tilt) if we knew what we're looking for...\nThe sentence has been modified to “…, rendering them unlikely to be observable only\nwith very high temporal resolution methods”.\nFigure 4: dotted orange line really hard to see...\nThe colour has been changed to brown to increase visibility.\nL209: add dot after \"rates\"\nAdded. Thanks.\nL233 - 238: should be observable in data - seen anywhere?\nThis comment has been addressed above.\nMany thanks for the comments.\nComments by Reviewer #2\nThis manuscript presents interesting results that highlight the role of thermos-elastic\ndeformation. I find of particular noteworthiness the results around the influence of thermal\nheterogeneity and how this produced corresponding heterogeneous viscoelastic deformation\nresponses. The work will be of significance to anyone in the field of volcano deformation.\nI see no major flaws in the method or analysis, and comments/suggestions for improvement\ncan be found in the attached file. These are minor and involve three main themes.\n1. A “toning” down of some of the claims in the introduction. This work included novel\naspects, but there have been several studies with comparable elements that have looked at\nnon-static thermally controlled deformation models of volcanic systems (a good few of which\nae cited by the authors).\nWe acknowledge that there are indeed several published models employing viscoelastic\nrheologies over timescales commonly used in volcano monitoring (e.g. Hickey et al.\n2016) as correctly pointed out by the reviewer. However, these studies consider the\ntemperature field and hence viscoelastic rheology to be static in time. We would\ntherefore like to emphasize that our study addresses a different aspect of non-static\ntemperature fields and therefore rheology by focusing on much longer timescales\nrelevant to magmatic system construction (~100-1000 ka), which we believe has not been\nexplored in the existing literature.\nTo address the reviewer's concern, we have rephrased the relevant section of the\nintroduction to better clarify the novelty of our study. The revised text now reads:\n“Given the temperature-dependence, the viscoelastic response of rocks subjected to\nstrain is largely governed by the long-term thermo-mechanical history of the crust39,69.\nWhile these timescales are much longer than those typically invoked in surface\ndeformation modelling, we anticipate that magmatic systems undergo different styles of\nground deformation at various stages of their lifecycle. This temporal evolution of\ncrustal rheology may be essential in explaining the great diversity of spatial and\ntemporal deformation behavior that is observed at different volcanoes worldwide40.\nYet, even state-of-the-art thermally controlled models of volcano deformation consider\nthe temperature field and hence viscoelastic rheology to be static, often relying on\nsimplified geometries of magma plumbing systems24,28. This is the first attempt to\nreconcile the vastly different timescales of ground movements and magmatic system\nevolution to explore these dependencies in space and time.” Lines 55-65\n2. I got there in the end, but I did find figure 4 initially difficult to understand, I may have not\nread thoroughly enough the opening few pages, or not realised the importance, but I think\nthis needs work to make sure thar the results are supported by the description of the\nmethodology, a really quick win would be referring back to Figure 1 when discussing Figure\nadded a clarification that refers back to Figure 1. The revised text now reads:\n“We model short-term deformation patterns over a time span of 10 years for\noverpressure sources at 5, 10 or 15 km depth and at different long-term evolutionary\nstages of magmatic systems (i.e. from 0 to 1 Ma in 200 ka increments; cf. Fig. 1).”\nAdditionally, the reference to Figure 1 has been included in the caption of Figure 4 to\nenhance clarity and support the description of the methodology.\n3. I have no issue with the underlying methodology, it is appropriate and understandable. I\nwould however like to see more commentary on the creation of the 800m x 200m pressurised\ncavity, or specifically more commentary on what it is not. I may be being oversensitive to\nthis, but the manuscript is presented as b the results being envisage as created from a distinct\nintrusion or equivalent, after a given period of magmatic activity. However, the result does\nnot really model the intrusion of anything, as the cavity is pre-formed in the model mesh and\nthen pressurised. Any comments therefore on the initial elastic response cannot be overly\ninterpreted, and the appropriateness of if representing an intrusion is also limited. The\nstrength, as mentioned previously is looking at the viscoelastic response, and the authors\nalready allude to this in lines 147-148. I would like to see it more explicitly expressed that the\ninitial elastic response does not represent any expected volcanic process, while emphasising\nthe observations on the thermoelastic response.\nWe appreciate the feedback and agree with the concerns raised. To address this point,\nwe have added a statement to the main text in Lines 193-199:\n“Although the relaxation behaviour in our model is governed by distribution of\nviscosities rather than single values, we do not model the injection process itself. The\nstep-change cavity pressurization does not reflect a specific igneous process and cannot\nbe considered equivalent to a forcing timescale due to magma injection. Instead, our\nstudy focuses on the time-dependent thermo-viscoelastic behavior of crustal rocks in\nresponse to a pre-formed, pressurized cavity.”\nand additional text in the methodology section (Lines 614-618):\n\"It is important to note, however, that we do not explicitly model the intrusion process\nthat leads to cavity formation. Consequently, the initial elastic response should not be\ninterpreted as resulting from a specific igneous process. Instead, our focus is solely on\nthe thermo-viscoelastic processes resulting from a pre-formed, overpressurized cavity. \"\nAdditionally, we have removed or revised any parts of the manuscript that implied the\ninitial elastic response was interpreted as an intrusion process.\nThese are minor points and mainly clarifications, I have made a few suggestions around\ncleaning up the text in the introduction, otherwise I enjoyed reading the manuscript which\ndetails an interesting and useful piece of work.\nThank you for the comments. We appreciate the constructive feedback.\nReviewer #3 (Remarks to the Author):\nThis manuscript takes a creative approach to the problem of reconciling different timescales\ninvolved in volcano deformation. The authors model 100s kyr thermal evolution in a crustal\nsection, which then provides initial conditions and defines static but spatially heterogeneous\nmaterial properties in a 1-10 year viscoelastic model of a single spheroidal intrusion. The\napproach merges notions of thermal conditioning and long-term intrusive flux with the\ncavity-type descriptions of magma chambers used to model active volcano deformation.\nThe numerical modeling is 2D, a mash-up of a commercial, non-open source FEM code\n(COMSOL) and an open source finite difference code. The authors identify some interesting\nand somewhat non-intuitive results that suggest a role for thermal maturity - and by\nextension the history of magmatic flux - in volcano deformation timeseries, including\nscenarios for which intrusion results in subsidence due to viscoelastic effects and cyclic\nbehavior. The authors apply these results in a qualitative way to data from volcanoes in the\nEast African Rift.\nI support publication – depending how the authors address the comment listed in the first\nbelow. But overall, it would be nice to see a more thorough parameter exploration to argue\nthe main points convincingly. The authors could also do a better job justifying various\nassumptions in the model, as they are not unique and future work will benefit from a clearer\nexposition.\nBefore I believe the results, I’d like to know whether the subsidence you observe in models is\nsimply a result of inconsistent model initial conditions.\nYou are modeling long-term assembly in a kinematic way, but the system has mechanical\nmemory and stresses associated with emplacement of intrusions have long-range influence.\nRucker et al., (2022) found that steady state cavity size varied with forcing period because\nviscoelastic creep is frequency dependent. This implies that the “size” of a magma chamber\ncan’t simply be imposed kinematically in time-dependent viscoelastic problems… it’s a model\nresult. My impression is that most viscoelastic chamber models do not address mechanical\ninitial conditions. But chamber volume impacts surface deformation, so this impacts the\nsignificance of your results. Similar issues relating to pre-stress have been long recognized in\nearthquake models (e.g., Day, 1982). We’re sortof behind the curve in volcanology.\nThank you for pointing out the findings of Rucker et al. (2022), which are certainly\nrelevant to the problem, and we now reference in Lines 75, 193, 470 and 610. Clearly,\nthe stress field and mechanical conditions are important factors to consider when\nmodelling magma injection location and reservoir size, but in many cases, these will be\noverprinted by crustal heterogeneities which are not possible to model. We agree that a\nlong-term goal should be a fully integrated model incorporating all these effects into a\nsingle modelling framework. However, the emphasis of this study is a new proof-of-\nconcept to explores a new link rather than the development of a new modelling\nframework. Thus, we treat these confounding factors illustratively in our modelling\napproach at present. We have added: ”In reality, the location and size of the source\nwould be controlled by the mechanical conditions and history-dependent stress\nfield63,68, but for illustrative purposes, we use an oblate source situated at depths of\neither 5, 10, or 15 km and test a range of source sizes.” to Lines 73-76.\nAs you correctly point out, reservoir size impacts surface deformation. We thus agree\nthat it is important to test whether our subsidence trends are impacted by reservoir size.\nIn the revised version of the manuscript, we present additional supplementary model\ncalculations with smaller and larger reservoir sizes (Fig. S7), showing that subsidence is\nalso observed in these scenarios, reinforcing our main finding.\nI’m also curious whether you’d see subsidence even in the absence of sill emplacement…\nsince you’re coupling two different codes (two entirely different numerical methods!), self-\nconsistency must be demonstrated. It could be that one needs to derive initial equilibrium\nmechanical conditions in an iterative manner (similar to what is done with gravitational\nloading from topography in elastic problems, e.g., Cianetti et al., 2012).\nAs stated in L581-585, our thermo-mechanical model implements a temperature field\ninto the benchmark approach of Hickey and Gottsmann (2014) that does not include\ngravity. Down-sagging of crustal rocks due to inconsistent mechanical equilibrium\nconditions is therefore not expected. Nevertheless, to test whether subsidence would\noccur without sill emplacement, we conducted simulations without applying any\nboundary load to the cavity. As expected, these simulations showed no displacement\nwhatsoever.\nOther comments:\nYou use Standard Linear Solid rheology but it’s a little unclear how you’re deriving a single\nrelaxation timescale (e.g., line 159).\nThe single value relaxation timescales are just example calculations to orient the reader\nwith regards to expected magnitudes of behaviour, bearing in mind the broad\nreadership of Nature Communications. Of course, the model domain is fully\nparameterized in term so relaxation timescales as a function of spatially variable\ntemperature. This is now stated more clearly in Lines 192-198.\nI am actually a bit skeptical of the need to use SLS rheology here despite several recent\npapers that advocate for it. SLS seems more a way to get sensible deformations from a\nphysically implausible model setup (constant step chamber pressure increase implies\nconstant magma flux into the chamber, does it not? Not even an oversimplified magmastatic\napproach to mass influx implies that…) rather than anything about crustal rheology.\nSpatially variable Arrhenius temperature dependence and time variable forcing in a Maxwell\nmodel will trade off with band-limited effects from the SLS model. Shouldn’t we use the\nsimplest model until real data dictate otherwise? Regardless, this all needs to be tightened\nup, explicit rheologic and mass balance assumptions stated etc.\nGiven that we are investigating general trends and pattern of surface deformation in an\nillustrative way, we would argue that the choice of Maxwell versus SLS rheological law\nis of secondary importance here. As you correctly point out, the question which\nrheology is more appropriate to real-world crustal rocks is currently discussed in the\nliterature. Differences between these approaches have been investigated in detail\nelsewhere (e.g. Head et al. 2021).\nNonetheless, we agree that some comparison of our deformation time series for different\nrheologies may be informative and added a supplementary figure (Fig. S5) comparing\nour results with SLS and Maxwell rheology. This analysis shows that the overall pattern\nand trends of surface deformation are conserved.\nWe now state explicitly the assumptions of constant step change injection scenario\n(Lines 192-198and Lines 614-618).\nAnd its not simply the relaxation time that governs deformation in viscoelastic systems as\ndescribed on line 159. It’s the ratio of relaxation time to forcing timescale (the Deborah\nnumber, for Maxwell rheology) that determines whether the surface deformation reflects\nelastic or viscous response (Rucker et al., 2022; Liao et al. 2023). I think you could do a\nmore convincing job with this problem if you developed the appropriate nondimensional\nnumbers and presented your results in this way. For example, it would be great to say\nsomething more concrete about how the spatial structure of the thermal field impacts the\nmechanical evolution (paragraph starting ln 223). What sets the lateral extent of the chamber\nor the thermal preconditioning here? It plays a big role. Then I’d seek to explicitly tie the\ninflux rate to the thermal preconditioning and then to the deformation. The latter is our main\nobservation, and it would be nice to know whether anything unique is encoded there. Could\nbe the answer is no.\nThank you for this stimulating discussion that lays out a compelling pathway for\nfurther analysis. A brief discussion of Deborah Numbers has been added to Lines 192-\n198 of the revised manuscript.\nWe can see that casting the results in terms of Deborah numbers could provide\nadditional insights. However, given that such calculations require knowledge of the\nforcing timescale, which in most modelled scenarios is treated as a constant cavity\nboundary load and therefore not representative of an injection timescale, this is beyond\nthe immediate scope of this work. However, thank you for pointing this out. It would\ncertainly be a good starting point for a follow-up study.\nIts worth noting that Liao et al (2023) has essentially resolved the dichotomy between older\n“shell” models and spatially variable thermoviscoelastic models of magma chamber\ndeformation. The “shell” can be precisely defined and derived analytically for each\nfrequency component of forcing, and one can associate the “shell” spatial extent with a\nDeborah number contour (even for step function pressure as modeled here). That’s a fairly\nbig result in my opinion.\nThis is an interesting finding. We cite this study in Lines 75, 193, 470 and 610.\nFinally, perhaps this is more of a philosophical comment, but I’m struck by some\nphenomenologic similarities between this manuscript and that in Karlstrom et al., (2010).\nFor example, Fig 4 in that paper doesn’t look so different to me than Fig 3 in the present\nmanuscript (appropriately blurring the eyes). Certainly the 2D stochastic approach to long-\nterm recharge and notions of thermal pre-conditioning as governing short-term dynamics are\nsimilar. There are so many benefits to the modeling framework developed here of course, but\nI suspect that the basic phenomenology of your results could be (and to some extent already\nhave been) more simply demonstrated.\nWe added a reference to your work (Line 54-56)\nGreat, provocative work regardless.\nThank you for the comments."}, {"round": 2, "reviewer_comments": {"reviewer_1": "The changes in response to my first round of comments are overall very satisfying. Thanks to the authors for thorough revisions!\nI think, however, that I lost a piece in my prior review transferring it to the online system. I think it's still worth including here. The manuscript might benefit from some discussion related to this, especially given the new text around lines 486-491.\nThe current claim that uniform uplift vs. uplift subsidence cycles are due to temperature / hot vs cold magmatic systems may be a bit strong. There have been other papers investigating various mechanisms unrelated to viscoelasticity, some of which already cited. More related to this work, though, Block et al. (2023), while using static rheology, suggested that complex spatial and temporal surface deformation pattern can be created over long time scales just by pressure variations over time (which are part of the model here too, at least in creating the subsurface structure) into a magma reservoir embedded in a viscoelastic region. They also demonstrate time dependent shifts in surface deformation, especially a phase lag between \"central\" and \"shoulder\" regions that they attribute to viscosity, and propose this as an alternative to the sombrero pattern observed at mid-crustal magma bodies. This is even for a relatively young (10kyrs?), but deeper system, the Socorro Magma Body, which seems to be at odds with the findings in this current paper. Perhaps that's due to the focus on the region of maximum uplift by Weber et al. I am not suggesting that the authors explore all temporal variations in their modeled signal - that's for later work that this paper opens up. But it seems that some more caution in the interpretation of their results is needed - clearly interpretations of surface deformation fields (often discontinuous and even missing regions with max deformation) will depend heavily on our prior knowledge of the system, and while we could interpret change in deformation over time as hot v. cold it could also be other reasons, such as viscous relaxation due to pressure reduction.\nOnce this aspect is addressed, I think this will make a fine contribution. I don't need to see the paper again.\nBlock, G.A., M. Roy, E. Graves, R. Grapenthin (2023), Pressurizing Magma Within Heterogeneous Crust: A Case Study at the Socorro Magma Body, New Mexico, USA, GRL, vol 50(20), e2023GL105689", "reviewer_2": "Thank you for taking the time to respond to my comments in a way which makes the manuscript, in my opinion, considerably easier to follow and understand.\nMy positive comments from the initial review still stand, and with the added clarity following revision, I think this is a very good piece of work and have no concerns with recommending publication in its revised form.", "reviewer_3": "I appreciate the considerate responses of the authors to reviewer comments, and feel the manuscript is much improved. Its a nice study.\nI'm still hung up on your justification of initial conditions, and subsequent discussion of the origin of subsidence. I don't think you addressed that concern completely.\nIts not simply that \"the location and size of the source would be controlled by the mechanical conditions and history-dependent stress field\" - the nature of the deformation response depends heavily on this. As you state, subsidence following cavity pressurization only occurs under certain influx conditions (which set the spatial distribution of temperature hence viscosity). Could you better predict when downflow will occur based on simple analysis of the viscous force balance? Given the neglect of gravity (which may be tenuous in the case of viscous flow... certainly justified for elastic response), seems like the existence of viscous response in downward direction could be straightforwardly predicted based on magnitude of overpressure versus viscosity structure.\nBut then you are assuming a uniform overpressure and zero initial stress. What about the magma injection that (presumably) is your preferred origin for intrusion pressurization? That must have itself flowed down pressure gradient into the chamber. Would that deeper pressure source balance the downward flow of hot host rock? Thats the type of \"pre-stress\" I was referring to. I think you need to more clearly contextualize the assumptions you're making about magma mass transport, even if thats not modelled explicitly. I'm guessing that this subsidence behavior is an end member at best."}, "author_response": "Reviewer #1 (Remarks to the Author):\nThe changes in response to my first round of comments are overall very satisfying. Thanks to the\nauthors for thorough revisions!\nI think, however, that I lost a piece in my prior review transferring it to the online system. I think it's\nstill worth including here. The manuscript might benefit from some discussion related to this,\nespecially given the new text around lines 486-491.\nThe current claim that uniform uplift vs. uplift subsidence cycles are due to temperature / hot vs cold\nmagmatic systems may be a bit strong. There have been other papers investigating various\nmechanisms unrelated to viscoelasticity, some of which already cited. More related to this work,\nthough, Block et al. (2023), while using static rheology, suggested that complex spatial and temporal\nsurface deformation pattern can be created over long time scales just by pressure variations over time\n(which are part of the model here too, at least in creating the subsurface structure) into a magma\nreservoir embedded in a viscoelastic region. They also demonstrate time dependent shifts in surface\ndeformation, especially a phase lag between \"central\" and \"shoulder\" regions that they attribute to\nviscosity, and propose this as an alternative to the sombrero pattern observed at mid-crustal magma\nbodies. This is even for a relatively young (10kyrs?), but deeper system, the Socorro Magma\nBody, which seems to be at odds with the findings in this current paper. Perhaps that's due to the\nfocus on the region of maximum uplift by Weber et al. I am not suggesting that the authors explore all\ntemporal variations in their modeled signal - that's for later work that this paper opens up. But it seems that\nsome more caution in the interpretation of their results is needed - clearly interpretations of surface\ndeformation fields (often discontinuous and even missing regions with max deformation) will depend\nheavily on our prior knowledge of the system, and while we could interpret change in deformation\nover time as hot v. cold it could also be other reasons, such as viscous relaxation due to pressure\nreduction.\nThank you for sharing this discussion. We now acknowledge the mechanism proposed by Block\net al. 2023 in Lines 44-45 of the revised manuscript and rephrased the text in the “Implications\nfor volcanology” section more carefully in Lines 487 and 499.\nOnce this aspect is addressed, I think this will make a fine contribution. I don't need to see the paper\nagain.\nBlock, G.A., M. Roy, E. Graves, R. Grapenthin (2023), Pressurizing Magma Within Heterogeneous\nCrust: A Case Study at the Socorro Magma Body, New Mexico, USA, GRL, vol 50(20),\ne2023GL105689\nThe Reference has been added.\nReviewer #2 (Remarks to the Author):\nThank you for taking the time to respond to my comments in a way which makes the manuscript, in my\nopinion, considerably easier to follow and understand.\nMy positive comments from the initial review still stand, and with the added clarity following revision,\nI think this is a very good piece of work and have no concerns with recommending publication in its\nrevised form.\nThank you.\nReviewer #3 (Remarks to the Author):\nI appreciate the considerate responses of the authors to reviewer comments, and feel the manuscript\nis much improved. Its a nice study.\nThank you.\nI'm still hung up on your justification of initial conditions, and subsequent discussion of the origin of\nsubsidence. I don't think you addressed that concern completely.\nIts not simply that \"the location and size of the source would be controlled by the mechanical\nconditions and history-dependent stress field\" - the nature of the deformation response depends\nheavily on this. As you state, subsidence following cavity pressurization only occurs under certain\ninflux conditions (which set the spatial distribution of temperature hence viscosity). Could you better\npredict when downflow will occur based on simple analysis of the viscous force balance? Given the\nneglect of gravity (which may be tenuous in the case of viscous flow... certainly justified for elastic\nresponse), seems like the existence of viscous response in downward direction could be\nstraightforwardly predicted based on magnitude of overpressure versus viscosity structure.\nAs you correctly point out, the occurrence of subsidence can indeed be predicted based on the\nviscosity structure and the magnitude of overpressure. We have now incorporated an additional\nanalysis to address this point. Specifically, we provide a new figure 8, which illustrates how the\ndifference in relaxation timescales above and below the pressure source, calculated as η/G for\nthe Standard Linear Solid (SLS) rheology, can be used to predict the conditions under which\nsubsidence will occur. When the contrast in relaxation times is large, subsidence is likely to\noccur, whereas uplift tends to dominate when the difference is smaller. As you mention above,\nthis could also be expressed in terms of overpressure, with the relaxation timescale given by η/σ,\nyielding a similar result.\nWe added the following discussion:\n“The occurrence of either subsidence or uplift can be predicted by analysing the difference in\nrelaxation timescales of crustal rocks, which is primarily governed by their viscosity structure.\nWhen there is a significant difference in relaxation timescales between the rocks above and\nbelow the pressure source, subsidence is likely, as the upper layer behaves elastically while\nviscous relaxation occurs in the lower layer (see Fig. 8). In contrast, when the difference is\nminimal, both layers undergo time-dependent deformation, resulting in net uplift. Thus, the\nrelative difference in relaxation timescales, rather than absolute values, is the key factor\ncontrolling the deformation behaviour.”\nThank you for this comment. We believe it strengthened the manuscript.\nBut then you are assuming a uniform overpressure and zero initial stress. What about the magma\ninjection that (presumably) is your preferred origin for intrusion pressurization? That must have itself\nflowed down pressure gradient into the chamber. Would that deeper pressure source balance the\ndownward flow of hot host rock? Thats the type of \"pre-stress\" I was referring to. I think you need to\nmore clearly contextualize the assumptions you're making about magma mass transport, even if thats\nnot modelled explicitly. I'm guessing that this subsidence behavior is an end member at best.\nThank you for the clarification. To address your concern, we have significantly expanded the\ndiscussion surrounding the modelling assumptions. You raise a valid point that the context of\nthese assumptions requires further elaboration, although this discussion must remain somewhat\nspeculative, as fully coupled magma injection models incorporating viscoelastic rheology are\nstill beyond the current state-of-the-art. The following paragraph has been added:\n“In this study, we do not assume any specific mechanism for pressure generation within the\nmodelled cavity. However, two widely accepted processes for magma pressurization are the\ninjection of fresh magma into a preexisting reservoir via a feeder dike and the accumulation of\nan exsolved volatile phase within the magma body72. For volatile-driven mechanisms, where\npressure is generated internally, we do not expect these processes to directly influence the\nrelaxation behaviour of the crustal rocks. However, magma transport through a feeder dike\nsystem could affect the stress field beneath the intrusion, potentially counteracting the\ndownward-directed flow of hot host rocks. More complex models are needed to explore the\ninteraction between magma injection and viscoelastic crustal rheology to determine whether\nthis hypothetical effect would be significant. Nonetheless, petrologically constrained magma\nrecharge timescales, typically ranging from days to months72, are shorter than the 10-year\nobservation window modelled in this study. Therefore, our results are applicable to scenarios in\nwhich the injection process has terminated.”\nThe comments are very much appreciated. Thank you."}, {"round": 3, "reviewer_comments": {"reviewer_3": "The authors have addressed my concerns sufficiently and I'm happy to see this published."}, "author_response": ""}]

Decision letter