Exploring the Impact of Thermally Controlled Crustal Viscosity on Volcanic Ground Deformation

Volcanoes undergoing unrest often produce displacements at the ground surface, providing an important window to interpret the dynamics of the underlying magmatic system. The thermomechanical properties of the surrounding host rock are expected to be highly heterogeneous, with key physical parameters having a strong dependence on temperature. Deformation models that incorporate nonelastic rheological behaviors are therefore heavily reliant on the assumed thermal conditions, and so it is critical to understand how the thermomechanical crustal structure affects the observed deformation field. Here, we use a series of thermo‐viscoelastic Finite Element models to explore how variations in thermal constraints (i.e., reservoir temperature and background geothermal gradient) affect surface displacement patterns when using the Maxwell and Standard Linear Solid (SLS) viscoelastic configurations. Our results demonstrate a strong variability in the viscoelastic deformation response when changing the imposed thermal constraints, caused by the partitioning of deformation and the dissipation of induced stresses. When using the SLS rheology, we identify that cumulative long‐term displacements can vary by over 20%, relative to a reference model with a reservoir temperature of 900°C and background geothermal gradient of 30 K km−1. The relative change increases to a maximum of 35% when thermal weakening of the Young's modulus is also considered. Contrastingly, the deformation patterns of the Maxwell rheology are governed by unbounded displacements and complete stress relaxation. Ultimately, we outline that uncertainties in the thermal constraints can have a significant impact on best‐fit source parameters (e.g., size and depth) and overpressure/volume‐change loading histories inferred from thermo‐viscoelastic models.

When modeling ground deformation patterns in volcanic regions, the assumption of purely elastic behavior presents a simple foundation for interpreting first-order characteristics of the deforming reservoir. However, in nature, pure deformation of a single type (e.g., elastic) does not exist; instead, different types of deformation occur simultaneously and with differing proportions (Burov, 2011). Whilst elastic models are capable of reproducing observed uplift and subsidence patterns, providing estimates of the source location, they often rely on pressure changes with unrealistic amplitudes (Del Negro et al., 2009;Masterlark et al., 2010;Newman et al., 2006). Inferred pressure changes can be ill-determined unless the source dimensions are independently constrained; a ubiquitous challenge for geodetic modeling. The presence of shallow or long-lived magmatic systems in volcanic regions can significantly perturb the geothermal gradient (e.g., Annen, 2011;Gelman et al., 2013;Karakas et al., 2017). Elevated temperatures in the vicinity of magmatic systems are expected to characterize the rheology of the middle and upper crust (e.g., de Silva & Gregg, 2014) and invalidate the elastic approximation by facilitating the nonelastic behavior of crustal materials. Elastic behavior is generally limited to the deformation of crustal materials at temperatures lower than those expected at the brittle-ductile transition, and is dependent on the strain-rate and crustal composition (Del Negro et al., 2009;Newman et al., 2001;Ranalli, 1995). In these regions, a component of viscous behavior is expected to heavily influence the observed spatiotemporal deformation field (Del Negro et al., 2009;Jellinek & DePaolo, 2003). Viscoelasticity combines an instantaneous elastic response with time-dependent viscous behavior, and is often incorporated into deformation models to account for the heterogeneous thermomechanical properties in the vicinity of a magmatic system, known to influence and partition the resulting deformation field (e.g., Del Negro et al., 2009;Gottsmann & Odbert, 2014;Gregg et al., 2012;Hickey et al., 2015Hickey et al., , 2016. Different methods can be used to incorporate a viscoelastic rheology into models of volcanic ground deformation, for example, by employing: (a) a viscoelastic half-space (e.g., Bonafede et al., 1986), (b) viscoelastic shells embedded within an elastic domain (e.g., Delgado et al., 2018;Dragoni & Magnanensi, 1989;Galgana et al., 2014;Newman et al., 2001;Segall, 2019), (c) two-layered media composed of an elastic layer underlain by a viscoelastic domain (e.g., Hickey et al., 2013;Pearse & Fialko, 2010;Sigmundsson et al., 2020;Yamasaki et al., 2018Yamasaki et al., , 2020, or (d) a thermo-viscoelastic regime (e.g., Gottsmann et al., 2020;Hickey et al., 2016;Morales Rivera et al., 2019).
Aside from thermo-viscoelasticity, many of these models require the viscosity of the viscoelastic domain to be defined a priori. Alternatively, the viscosity of a viscoelastic shell can be determined by assuming a source loading function and solving for the rheology, based on the goodness-of-fit of the model results (e.g., Delgado et al., 2018;Novoa et al., 2019). Viscosity values in the range of 10 16 -10 19 Pa s are often attributed to viscoelastic shells (e.g., Dragoni & Magnanensi, 1989;Jellinek & DePaolo, 2003;Newman et al., 2001) to represent the local ductile wall-rock surrounding the magmatic system. However, for thermo-viscoelastic studies, the imposed thermal regime (i.e., reservoir temperature and background geothermal gradient) depends on both the volcanic center and the availability of data constraints. Whilst the reservoir temperature can be reasonably estimated from petrological studies of eruptive products, a geothermal gradient of 30 K km −1 is often assumed. However, several studies also consider elevated or depth-dependent geothermal gradients, as shown in Table 1.
Viscoelastic media exhibit distinct time-dependent behaviors determined by stress and strain conditions. As a result, viscoelastic media produce different spatiotemporal deformation patterns, dependent on whether the source deformation processes are represented by a stress-controlled pressure change (ΔP) or a strain-controlled volume change (ΔV); two modes of deformation which can be considered equivalent when using an elastic rheology (Head et al., 2019). With a variety of modeling approaches and rheological models available, it is important to ensure that the deformation model is an accurate approximation of reality. If an inappropriate rheology is employed, the inferences of the underlying deformation process may be implausible, despite the construction of models that can satisfy goodness-of-fit criteria. We extend the work of Head et al. (2019)  and background geothermal gradients (ranging from 20 to 50 K km −1 ), and investigate the effect of the resultant temperature-dependent viscosity distributions on the modeled ground displacements. By imposing assumed reservoir dynamics as boundary conditions to the rheological response, we present generalized results that can be applied to different volcanic settings to understand how uncertainties in imposed thermal constraints can affect surface deformation in thermo-viscoelastic models.

Numerical Modeling
We use COMSOL Multiphysics® (v5.3a) to construct and solve finite element models that couple the structural mechanics and heat transfer modules, and evaluate the spatiotemporal deformation patterns relating to thermo-viscoelastic regimes. By applying thermal constraints, the resultant temperature field, composed of a background geothermal gradient and the thermal perturbation owing to the modeled reservoir, is used to determine a temperature-dependent viscosity distribution. In this investigation, we compare and contrast the temperature-and time-dependent behaviors of linear viscoelastic models, determined by the stress and strain conditions produced by each deformation mode (Figure 1a), and include the commonly used Maxwell and SLS configurations (Figures 1b and 1c) in our analysis.

Model Setup and Geometry
The model space is adapted from the benchmarked approach outlined in Hickey and Gottsmann (2014), shown in Figure 2, with r-and z-dimensions of 30 km. The deforming magma reservoir is represented by a finite spheroidal cavity, with a radius of 1,500 m and centered at a depth of 5 km. In parallel model setups, the reservoir wall is allocated either an overpressure (ΔP) or a prescribed displacement (Δ⍺) deformation mode; both are applied normal to the reservoir wall. To ensure agreement between the overpressure (ΔP) and volume-change (ΔV) models, we use elastic models to derive a prescribed displacement (Δ⍺) that produces the same vertical displacement at the ground surface (z = 0) as the overpressure models (ΔP), simulating the process of geodetic inversions. It should be noted that the volume-change models require a greater displacement on the reservoir boundary than is produced by the overpressure models; an effect imparted by the free surface of the finite model domain (i.e., not a full-space; Figure 2a). For simplicity, we impose a uniform displacement on the reservoir walls (producing an isotropically expanding source) as a first-order approach for implementing a strain-controlled volume-change, a reasonable assumption for a spherical deformation source. Ensuring a complete match between the overpressure and volume-change models, in terms of induced stresses and displacements on the reservoir wall, is out-of-scope for this study. We use a fixed overpressure of 10 MPa whilst the prescribed displacement (and therefore volume-change) HEAD ET AL.  varies on a model-by-model basis, as they are dependent on the imposed thermal constraints. Vertical displacement time series are evaluated at the point on the free surface that lies directly above the center of the modeled reservoir (i.e., r = z = 0).
Traditionally, with elastic modeling, the inferred temporal reservoir evolution, and hence the associated overpressure or expansion requirements, directly reflects the profile of a deformation time series. We investigate the deformation response associated with a constant loading function to highlight the viscoelastic effect. A constant (time-invariant) overpressure or volume-change may be used to represent a reservoir that has reached equilibrium following a prior intrusion, or a scenario where the influx and outflux of volatiles or magma, or combination thereof, are equal. In a later section, we discuss the implications for a ramped loading function, commonly used in volcano deformation studies, which may represent a continued injection or accumulation of magma, exsolution of volatiles, or combination thereof. The influence of the temperature-dependent viscosity structure can be extended to other loading functions (such as exponential, logarithmic, or trapezoidal), due to the Boltzmann superposition principle. We model the deformation over a time period of 10 years with a temporal resolution of at least 0.01 year. This timescale is directly relevant to volcano monitoring timeframes (Phillipson et al., 2013), and is of sufficient length to observe the time-dependent behavior of rocks with a viscosity of less than 10 19 Pa s. Our models do not explicitly account for the effect of gravitational loading, as we investigate the deformation produced when loading away from an equilibrium state. We implement static deformation modes in this investigation and so neglect any second-order dynamical pressure-volume relationships, such as the changes in overpressure induced by a change in volume (e.g., Gregg et al., 2013).

Temperature Distribution and Viscosity
The temperature distribution is computed numerically using the steadystate heat conduction equation: where k is the thermal conductivity, T is the temperature, and Q is the crustal volumetric heat production, which is assumed to be zero. In this study, the thermal conductivity and specific heat capacity are constant. Any uncertainties in the values of the thermal conductivity and specific heat capacity are negligible for our calculation; they only affect the timescale of heat conduction (Hickey et al., 2015;Morales Rivera et al., 2019) and we assume that the system has reached thermal equilibrium. The temperature distribution is constrained by fixing the surface temperature, s T , and background geothermal gradient, dT dz , which is then perturbed locally by the reservoir temperature, r T . Herein, "geothermal gradient" refers to the background geothermal gradient imposed in each model variant, whilst "thermal regime" refers to any given combination of reservoir temperature and geothermal gradient. The thermal constraints are shown in Table 2. By ascribing a constant temperature to the source boundary, the reservoir walls are considered to be acting as a heat source (Del Negro et al., 2009).
We calculate the thermally controlled viscosity structure by using the Arrhenius formulation: HEAD ET AL.
10.1029/2020JB020724 4 of 22 The springs and dashpots are characterized by the shear modulus, G, and shear viscosity, η, respectively. In the SLS model, the shear modulus is split across the elastic and Maxwell arms by the fractional moduli, μ 0 and μ 1 . Full behavioral descriptions of the Maxwell and SLS models under different stress and strain conditions can be found in Christensen (1982), Crawford (1998), or Head et al., (2019. where D A is the Dorn parameter, 1 × 10 9 Pa s (Del Negro et al., 2009;Gregg et al., 2012;Hickey et al., 2016), A E is the activation energy, 1.3 × 10 5 kJ mol −1 (Meissner & Tanner, 1992;Ranalli, 1995), R is the universal gas constant, and T is the temperature. Through this relation, elevated temperatures in the vicinity of the modeled reservoir result in reduced viscosities ( Figure 2d). To ensure that excessively low viscosities are not produced at depth, the geothermal gradients are applied to a maximum temperature of 950°C, a conservative estimate for the deep crust in volcanic regions (e.g., Hickey et al., 2016;Schutt et al., 2018). For geothermal gradients of 40 and 50 K km −1 , this results in an isothermal and isoviscous region, of 10 14.4 Pa s, at the base of the model. Without this limit, geothermal gradients of 40 and 50 K km −1 could produce viscosities as low as 10 12.2 and 10 11.2 Pa s, respectively.

Thermomechanical Heterogeneity
We introduce mechanical heterogeneity within the models by adjusting the Young's modulus, E, in the vicinity of the modeled reservoir to account for thermal weakening; where the Young's modulus can decrease HEAD ET AL.
10.1029/2020JB020724 5 of 22 Figure 2. Diagrams of the model setup. Geometry, boundary conditions (a), and thermal constraints (b) of the 2D-axisymmetric finite element model, with a spherical cavity representing the deformation source. To prevent excessive temperatures at depth, the imposed geothermal gradients are conservatively limited to a maximum temperature of 950°C (e.g., Schutt et al., 2018 by over 50% at temperatures above 900°C (Bakker et al., 2016;Smith et al., 2009). Following the approach of Zhan et al., (2019), we use an exponential smoothing function of the form: where T is the steady-state conductive temperature distribution determined from the imposed source temperature, r T , and geothermal gradient, dT dz . Away from the modeled reservoir (∼10 km), this relationship returns an unaltered Young's modulus, E (Figure 2e).
The viscous component of deformation is commonly assumed to be incompressible; such that volumetric strains are considered to be purely elastic and viscoelastic deformation is represented in terms of the deviatoric strain components (Segall, 2010). This allows the bulk modulus, K, to behave elastically whilst the HEAD ET AL.  shear modulus, G, behaves in a viscoelastic manner (Del Negro et al., 2009;Folch et al., 2000). The bulk and shear moduli are determined from the Young's modulus used in each model, using standard formulation.

Results
Here we compare the resultant vertical displacement time series for the contrasting deformation modes (i.e., ΔP or ΔV) within the Maxwell and SLS thermo-viscoelastic models. The displacements are evaluated at the point on the ground surface that directly overlies the centroid of the modeled reservoir (i.e.,   0 r z ).

Sensitivity to Thermomechanical Heterogeneity
Prior to investigating the effect of thermo-viscoelasticity on ground deformation, we consider the sensitivity of the deformation models (overpressure, ΔP, and volume-change, ΔV) to each temperature-dependent component. The degree of thermomechanical heterogeneity, and therefore model complexity, is increased in a stepwise manner, by considering: The initial model configuration consists of a spatially uniform Young's modulus (20, 30, and 40 GPa) and prescribed viscosity (10 16 , 10 17 , and 10 18 Pa s), in a similar approach to Head et al. (2019). Through Equation 2, the viscosity values used in these isoviscous models relate to temperatures of ∼700, ∼575, and ∼480°C, respectively. This initial configuration is then perturbed using a reference thermal model, with a reservoir temperature of 900°C and a geothermal gradient of 30 K km −1 . We examine how a temperature-dependent Young's modulus and temperature-dependent viscosity distribution affect deformation in separate, but complementary, steps before combining both components in a final stage. Resultant time series for the sensitivity analysis of the overpressure (ΔP) and the volume-change (ΔV) models can be seen in Figures

Overpressure
With the overpressure (ΔP) deformation mode, the reservoir boundary experiences a constant force directed away from the center of the reservoir (Figure 1a). In response to the applied stress, the viscoelastic media undergo a time-dependent increase in strain, known as creep (e.g., Christensen, 1982;Crawford, 1998;Ranalli, 1995), which is characterized differently for the Maxwell and SLS models.
For the Maxwell viscoelastic configuration, the isoviscous (Stages 1 and 2a) models demonstrate linear, unbounded displacements in response to the imposed overpressure (e.g., Head et al., 2019), with the gradient determined by the viscosity of the model domain ( Figure 3a). However, incorporating a temperature-dependent viscosity structure (Stages 2b and 3) enables the rheology to produce a rate-decreasing deformation response with time, the magnitude of which is dependent on the Young's modulus ( Figure 3a). By normalizing the time series, shown in Figure 3b, the displacements corresponding to the Stage 3 models are ∼6 times greater than the elastic response by the end of the 10-year modeling period, but are still increasing. Including a temperature-dependent Young's modulus (Stages 2a and 3), decreasing the strength of the crustal rock in the vicinity of the modeled reservoir, makes negligible difference to the absolute magnitude of the viscous displacement ( Figure 3a) but the displacement relative to the corresponding elastic model is increased (Figure 3b).
In response to the imposed overpressure, the SLS rheology produces rate-decreasing (asymptotic) deformation (e.g., Head et al., 2019). The isoviscous models (Stages 1 and 2a) demonstrate that the magnitude of HEAD ET AL.  viscous displacement is independent of the viscosity, with the 10 16 and 10 17 Pa s models converging to the same displacement value over the modeled period ( Figure 3c). Normalizing the time series (Figure 3d) identifies an asymptote of ∼1.71 times that of the elastic response for the models with a uniform Young's modulus (e.g., Del Negro et al., 2009;Head et al., 2019). Incorporating a temperature-dependent Young's modulus (Stage 2a) increases the magnitude of deformation observed, increasing the asymptote to a factor of ∼1.79 (Figure 3d). Using a temperature-dependent viscosity structure (Stages 2b and 3; Figure 2d) increases the initial rate of displacement, due to the range of viscosities produced, but also decreases the magnitude of viscous deformation with respect to the isoviscous models (Figures 3c and 3d).

Volume-Change
The volume-change (ΔV) deformation mode applies a uniform displacement directed away from the center of the reservoir. In response to the strain condition the viscoelastic media experience a time-dependent dissipation of stress, known as relaxation (e.g., Christensen, 1982;Crawford, 1998;Ranalli, 1995), which is characterized differently for the Maxwell and SLS viscoelastic configurations.
Unlike the creep behavior in response to overpressure, the relaxation behavior of the Maxwell and SLS viscoelastic configurations both exhibit an exponentially decreasing deformation within the isoviscous models (Stages 1 and 2a), as shown in Figures 4a and 4c, with an asymptote that is clearly defined by the 10 16 and 10 17 Pa s time series. The absolute magnitude of the resulting displacement is dependent on the Young's modulus, whilst the rate of decay is dependent on both the Young's modulus and the viscosity of the model domain (e.g., Head et al., 2019). However, the Maxwell and SLS configurations have different relaxation times, which control the decay rate of the displacement, and deformation asymptotes. For a uniform Young's modulus, the Maxwell rheology displays an asymptote of ∼0.65 times that of the elastic response (Figure 4b), whilst the SLS rheology converges to a value of ∼0.84 (Figure 4d). Incorporating a temperature-dependent Young's modulus (Stage 2a) has a minor impact on the absolute magnitude of the relaxation response for both viscoelastic configurations, as seen in Figures 4a and 4c, however, the displacement asymptotes are increased to a factor of ∼0.77 and ∼0.89 for the Maxwell and SLS rheologies (Figures 4b and 4d), respectively.
The inclusion of a temperature-dependent viscosity distribution (Stages 2b and 3) has a large effect on the displacement time series. The Maxwell rheology ( Figure 4a) exhibits an initial increase in displacement prior to a long-term decay signal, which arises due to the presence of low viscosities in the vicinity of the reservoir and the corresponding timescales for stress relaxation. In contrast, the SLS configuration produces a short-term, low-magnitude decrease in displacement followed by a partial recovery ( Figure 4c). Adding temperature-dependence to the Young's modulus (Stage 3) results in a reduction of the timescales relative to the Stage 2a time series (e.g., Figure 4b). The SLS rheology produces displacements with reduced magnitude and appears to converge to a value of ∼0.98 times to that of the elastic response, with a minor difference for the Stage 2a model. In contrast, the modeling period is not sufficiently long to identify the limit of decay for the Maxwell configuration.

Viscous Timescales
The timescales of the creep and relaxation behaviors displayed for both viscoelastic models are controlled by a characteristic timescale, which is given by the ratio of the viscosity, , to the shear modulus, G (e.g., Christensen, 1982;Head et al., 2019;Ranalli, 1995): In response to the imposed overpressure, the Maxwell creep behavior is linearly dependent on  , whereas for the relaxation behavior the dependence is exponential. The creep and relaxation behavior for the SLS rheology is also exponentially dependent on  , as well as an additional dependence on the inverse product of the fractional moduli (i.e.,   and 4 are determined by both the viscosity and the elastic moduli, however, the elastic moduli play a secondary role due to the relative magnitude of variation; the low-and high-Young's modulus are of similar magnitude, whereas there are two orders of magnitude between the low-viscocity and high-viscosity values. For example, the viscosity values considered for the isoviscous models (10 16 -10 18 Pa s, relating to temperatures between ∼700 and ∼480°C) produce deformation responses on an ∼annual-scale to decadal-scale, converging to a constant value after sufficient time has elapsed (Figures 3 and 4).
Viscous effects are usually assumed to be limited to long-term deformation responses, and negligible over the short-term, which may be the case for models that utilize a model space viscosity of greater than 10 17 Pa s. However, the variation in crustal viscosity produced by the thermomechanical models allows viscoelasticity to capture both the rapid deformation response of the hotter, more-ductile crustal material in the vicinity of the magmatic reservoir, as well as the long-term response of the cooler material near the ground surface. Approaching the ground surface, where the viscosity of the crustal rock exceeds 10 20 Pa s (e.g., Figure 2d), the deformation timescales become sufficiently long (∼10 3 years) to approximate the elastic response over short monitoring periods. However, when considering a system with a heterogeneous crustal viscosity structure, such as that which exists in volcanic regions, it becomes very difficult to disentangle the rheological timescales and response from a time-varying deformation source; a persistent challenge faced by volcanic ground deformation models.

Variations in Thermal Regime
Results of the sensitivity tests demonstrate that incorporating a temperature-dependent viscosity distribution has a major impact on the predicted ground deformation timeseries. In this section, we explore how varying the reservoir temperature and geothermal gradient, which in turn alters the viscosity distribution, influences the resultant surface displacements. We explore a suite of thermomechanical configurations; with reservoir temperatures that reflect felsic, intermediate, and mafic magma compositions (700, 900, 1100°C, respectively) and a range of geothermal gradients (20-50 K km −1 ). Using Equation 2, the reservoir temperatures of 700, 900, and 1100°C produce minimum viscosity values of 10 16 , 10 14.6 , and 10 13.9 Pa s, respectively, at the interface with the host rock ( Figure 5). A background Young's modulus of 30 GPa decreases with temperature toward the reservoir, by ∼36-56% (Equation 3; Figure 2e) depending on the imposed thermal constraints, with high reservoir temperatures and low geothermal gradients resulting in the highest reductions.
HEAD ET AL.
10.1029/2020JB020724 11 of 22 Figure 5. Changes in temperature and viscosity with depth along the central axis (  0 r ) of the model. Three sets of lines are displayed, corresponding to reservoir temperatures of 700, 900, and 1100°C. The geothermal gradient is a primary control on the ambient host rock viscosity, varying by over three orders of magnitude beneath the modeled reservoir for a given reservoir temperature (e.g., 900°C, shown by a bold line of each color). In contrast, the temperature and viscosity distribution above the reservoir is primarily controlled by the reservoir temperature. The gray shaded region represents the source cavity.

Overpressure
The overpressure (ΔP) deformation mode applies a stress to the reservoir wall, with the resultant deformation, on both the reservoir and at the ground surface, determined by the viscoelastic rheology used. The strain incurred by the Maxwell configuration is unbounded and linearly related to the stress imposed, and therefore results in continued deformation of the reservoir whilst the overpressure is maintained. This behavior is characteristic of a viscoelastic fluid (Banks et al., 2011;Marques & Creus, 2012). In contrast, the strain response of the SLS rheology is limited, related to the weightings of the elastic and viscoelastic arms (  0 and  1), and therefore the displacements on both the reservoir and at the ground surface are restricted.
The relative magnitude and rate of the surface displacements produced by the viscoelastic models (Figure 6) are governed by the viscosity distribution of each thermal regime ( Figure 5).
In response to the imposed overpressure, the magnitude of the resultant displacement for the Maxwell rheology increases with reservoir temperature, as shown in Figure 6a, however, the geothermal gradient plays a primary role in determining the nature of the time series. Models with a low geothermal gradient (20 and 30 K km −1 ) appear to produce stable, rate-decreasing uplift, whilst higher geothermal gradients (40 and 50 K km −1 ) result in long-term post-uplift subsidence. In contrast, the SLS rheology produces a consistent deformation response across the range of thermal regimes considered, displaying rate-decreasing uplift. The magnitude of this displacement increases by both increasing the reservoir temperature and decreasing geothermal gradient (Figure 6b). In isoviscous media, the SLS rheology exhibits a fixed deformation asymptote with respect to the corresponding elastic response (∼1.79, Figure 3d), but this deformation limit appears to be strongly controlled by the heterogeneous viscosity structure. In Figure 6b, the apparent asymptote for a geothermal gradient of 40 K km −1 increases from a factor of ∼1.25-∼1.4 as the reservoir temperature is increased from 700 to 1100°C, whereas these limits are exceeded by models with lower geothermal gradients. HEAD ET AL.  The rate of the creep behavior exhibited by both viscoelastic configurations is dependent on the viscosity of the surrounding host rock ( Figure 5). Low viscosities in the vicinity of the reservoir enable high rates of deformation which, following the application of the overpressure, induces a component of near-instantaneous deformation for both viscoelastic rheologies. As an increased reservoir temperature decreases the minimum viscosity on the reservoir wall, reducing the characteristic timescale (i.e., Equation 4); this results in the near-instantaneous deformation component to increase in magnitude with increased reservoir temperatures ( Figure 6). Furthermore, as the strain response of the Maxwell rheology is unbounded, unlike the SLS configuration, the magnitude of the displacement increases as the reservoir temperature is increased, due to the decrease in viscosity near the reservoir.

Volume-Change
For the volume-change (ΔV) deformation mode, the reservoir wall undergoes an instantaneous strain at t = 0, which remains fixed for the duration of the modeling period such that the reservoir itself undergoes no further deformation. As a result, the surface displacements observed in the viscoelastic models (Figure 7) arise purely from the time-dependent dissipation of the induced stress field, the rate of which is controlled by the viscosity distribution of each thermal regime. Stresses within a Maxwell viscoelastic medium are capable of relaxing fully, a behavior that is characteristic of a fluid (Banks et al., 2011;Marques & Creus, 2012), whereas the dissipation is limited for the SLS rheology. As the rate of the relaxation response for both viscoelastic configurations is dependent on the viscosity of the surrounding host rock, using an increased reservoir temperature, which decreases the minimum viscosity on the reservoir walls, has a large effect on the resultant deformation.
The influence of the thermal regime is clearly seen in the displacement time series of the Maxwell rheology (Figure 7a), where the low reservoir temperature (700°C) produces a long-term uplift signal whilst the higher temperatures (900 and 1100°C) result in long-term subsidence. This temporal trend reflects the degree of stress relaxation that has occurred, where the timescales associated with the 1100°C reservoir temperature HEAD ET AL.  are sufficiently short to allow for significant dissipation. The geothermal gradient has a secondary control on the deformation response. For a given reservoir temperature, varying the geothermal gradient changes the magnitude of the deformation response; low geothermal gradients produce the greatest displacements. The variation in magnitude is related to changes to the ambient viscosity distribution ( Figure 5), with higher geothermal gradients reducing the characteristic timescales, and therefore limiting the magnitude of deformation attained. A similar effect is also seen for the SLS rheology, as shown in Figure 7b; whereby the reservoir temperature has the greatest impact on the temporal trend of the resultant time series, and the geothermal gradient affects the magnitude of the displacement. The majority of the considered thermal regimes result in long-term displacements that are reduced with respect to the elastic response, except when coupling a low reservoir temperature (700°C) with lower geothermal gradients (20 and 30 K km −1 ). Whilst the deformation response of the SLS rheology varies with changes to the thermal regime, the magnitude of the resultant long-term displacements differ from the corresponding elastic response by less than 5%. Figure 7 demonstrates that the time series produced by the 1100°C reservoir temperature are "accelerated" versions of those produced by the lower reservoir temperatures, due to the reduced characteristic timescales associated with the decreased viscosities. This acceleration is accompanied by a reduction in displacement. As a result, this suggests that if the modeling period were sufficiently long, the lower reservoir temperatures would produce time series with the same temporal pattern displayed by the 1100°C reservoir temperature, but with increased magnitude.

Discussion
Our results demonstrate the differences in deformation patterns produced by the overpressure (ΔP) and volume-change (ΔV) deformation modes for each viscoelastic rheology, Maxwell and SLS, as highlighted by Head et al., (2019). Further to this, we illustrate the impact of the thermal regime, and therefore heterogeneous viscosity structure, on the displacements produced. In this section, we analyze the resultant deformation patterns of the thermomechanical viscoelastic models, identifying the drivers behind the observed deviation from the viscoelastic behaviors usually presented, and comment on the implications for modeling ground deformation.

Thermomechanical Strain Partitioning
Whilst the reservoir temperature used, and therefore the minimum viscosity of the surrounding host rock, is a major control on the deformation timeseries observed in Figures 6 and 7, it can be seen that the geothermal gradient affects the relative magnitude of the resultant displacements. This arises due to the temperature, and viscosity, gradients induced by the different thermal regimes (e.g., Figure 5), which play an important role in the partitioning of deformation. As the volume-change (ΔV) deformation mode emplaces a fixed strain to the reservoir wall, it does not deform any further. Therefore, differences in the viscosity gradients proximal to the modeled reservoir affect the timescales associated with the relaxation of the induced stress field, and so increased geothermal gradients enable a greater degree of stress dissipation to occur. In contrast, the overpressure (ΔP) deformation mode allows the reservoir wall to deform freely in response to the imposed overpressure, dependent on the characteristics of the rheology used. Foremost, Figure 6 demonstrates that increasing the geothermal gradient results in a change from rate-decreasing uplift to subsidence for the Maxwell rheology, whereas changing the geothermal gradient affects the magnitude of displacement for the SLS configuration. Using a reservoir temperature of 900°C as an illustrative example, Figure 8 displays the total displacement along the reservoir wall at t = 10 year, normalized by the elastic displacement at the top of the reservoir (i.e., z = −3.5 km), for the Maxwell and SLS rheologies (Figure 8a). The plots demonstrate the change in strain partitioning as the geothermal gradient increases. Increased temperatures and reduced viscosities at depth, associated with the higher geothermal gradients (40 and 50 K km −1 ), allow the Maxwell rheology to partition elevated displacements beneath the reservoir (Figure 8a) where the ambient viscosity is low ( Figure 5). This is a consequence of the unbounded strain response of the Maxwell model. Even for the lowest geothermal gradient, 20 K km −1 , the deformation at the base of the reservoir is over 10 times greater than the elastic displacement. Figure 8b illustrates the displacements at the ground surface for geothermal gradients of 20 and 50 K km −1 , evaluated through time. Notably, these profiles highlight both the magnitude of displacements produced by the Maxwell rheology, as well as the spatial footprint of the HEAD ET AL.

10.1029/2020JB020724
14 of 22 subsidence produced in the presence of a high geothermal gradient (Figure 6a), with respect to the elastic model (i.e., t = 0). The elevated displacements at the base of the reservoir (Figure 8a) induce a significant viscous flow which, due to the conservation of mass, requires the region above the modeled reservoir to subside. The surface displacements (Figure 8b) also show that the viscous flow affects horizontal displacements toward the edge of the model. As the pronounced subsidence arises due to a characteristic fluid behavior (i.e., unbounded strain response), we caution against the use of the Maxwell rheology.
Whilst the SLS configuration also displays increased partitioning toward the base of the reservoir (Figure 8a) with increasing geothermal gradients, the magnitude of the displacements is much lower compared to the Maxwell model. As the strain response of the SLS rheology is limited, viscous flow beneath the reservoir is negligible. The initial decrease in normalized displacement with depth reflects the free surface effect, which would continue to decrease toward the base of the reservoir in an elastic medium. This effect HEAD ET AL.

10.1029/2020JB020724
15 of 22 Figure 8. Displacement partitioning along the reservoir wall for the Maxwell and Standard Linear Solid thermoviscoelastic models using the overpressure (ΔP) deformation mode, and associated displacement profiles evaluated at the ground surface through time. (a) Normalized displacements along the reservoir wall at t = 10 year, using a reservoir temperature of 900°C. The plots are normalized to the elastic displacement at the top of the reservoir (i.e., z = −3.5 km), and illustrate the strain partitioning that occurs with changes to the geothermal gradient. (b) Normalized displacement transects across the ground surface for a reservoir temperature of 900°C, and geothermal gradients of 20 and 50 K km −1 (top and bottom row, respectively). The vertical (Uz) and horizontal (Ur) displacements are normalized to the corresponding maximum elastic displacement (i.e., at t = 0). These plots highlight the magnitude of the displacements produced by the Maxwell rheology, including the spatial footprint of the subsidence signal observed in Figure 6a.
is not clearly observed for the Maxwell rheology (Figure 8a) due to the magnitude of the displacements. As noted by Gottsmann and Odbert (2014) when modeling a deep-crustal "hot zone", an increasing component of reservoir deformation is accommodated by more-viscous material toward the base of the reservoir. As such, reduced displacements are observed at the ground surface (Figure 6b) as the geothermal gradient is increased. Based on this viscosity-induced partitioning, inferences of the best-fit source parameters (e.g., size and depth) and overpressure loading history from thermo-viscoelastic models may vary appreciably depending on the reservoir temperature and geothermal gradient used.

Application to an Increasing Overpressure Load
Whilst the results of this study use a constant reservoir loading, featuring an abrupt step-change at the onset (i.e., t = 0), volcanic ground deformation models are generally used to investigate increases in displacement, relating to periods of unrest, and often consider a time-dependent reservoir loading function (e.g., Cabaniss et al., 2020;Delgado et al., 2018;Gottsmann et al., 2020;Hickey et al., 2016;Le Mével et al., 2016). Therefore, it is important to place this study within a similar context. Based on the magnitude and variation of displacements produced by the Maxwell rheology (Figures 6a and 7a), we limit this analysis to the SLS configuration, which produces a consistent, and limited, deformation response with changes in thermal regime. Due to the small variation in displacements (less than ∼5%; Figure 7b) in response to the volume-change (ΔV) deformation mode, we consider only the overpressure (ΔP) deformation mode. An alternative analysis is considered for Maxwell viscoelasticity in the following section, based on a more common implementation.
We consider a linearly increasing overpressure load of 1 MPa yr −1 as an illustrative example, with the resultant displacements displayed in Figure 9a. Over the course of the modeling period, the viscoelastic timeseries become distinct from the elastic responses (dashed lines), with the differences in the displacements increasing with time. This is a result of the Boltzmann superposition principal, whereby each increment of load produces an independent and additive contribution to the total displacement. The initial rates of HEAD ET AL.  for each of the thermal regimes considered, relative to a reference reservoir temperature of 900°C and geothermal gradient of 30 K km −1 . The long-term displacements across the different thermal regimes can potentially vary by up to 35% with a temperature-dependent Young's modulus ( Td E ). The gray region illustrates the component of the variability that results from the viscosity distribution, by using a uniform Young's modulus (30 GPa), with a potential misfit of greater than 20%. displacement are nonlinear, reflecting the high rates of displacement produced by the low viscosity of the surrounding host rock (Figure 6b), but the magnitude of the displacements are low. After sufficient time has elapsed, the rate of the viscoelastic timeseries becomes linear, relating to the displacement limits observed in Figure 6b.
To identify the variability of the deformation response with respect to the imposed thermal regime (Figure 9b), we normalize the displacements relative to a reference model. In this case, we use the 900°C reservoir temperature with a geothermal gradient of 30 K km −1 . By transforming the timeseries in this way, Figure 9b demonstrates that greatest variability is produced when coupling a low reservoir temperature with a high geothermal gradient (700°C and 50 K km −1 ; decreased relative to the reference) or a high reservoir temperature with a low geothermal gradient (1100°C and 20 K km −1 ; increased relative to the reference). As a result, the long-term displacements can vary by up to 35% across the range of thermal regimes considered. For a given reservoir temperature, changing the geothermal gradient can impact the resultant displacements by ∼15%. If the rate of overpressure loading is sufficiently high, or is sustained over a prolonged period, then the variability becomes more apparent due to the magnitude of the displacements. A component of this variability arises from the temperature-dependent Young's modulus, as it affects the partitioning of deformation. Using a uniform Young's modulus therefore isolates the viscous component of the observed displacement variability, and is shown by the gray shaded regions in Figure 9b. This reduces the variability across the thermal regimes with the long-term displacements varying by over 20%, compared to the ∼35% variability with a temperature-dependent Young's modulus. For a given reservoir temperature, changing the geothermal gradient results can vary displacements by ∼10%. The resultant viscosity distribution therefore plays a key role in the potential variability of the ground displacements produced.
Given that ground deformation models seek to constrain the characteristics and temporal evolution of an underlying source, uncertainties in model parameters are translated to the inferences of the deformation source. This is demonstrated by the displacement variability shown in Figure 9b, where uncertainties in the thermal parameters (i.e., reservoir temperature and geothermal gradient) will directly affect the inferred best-fit source properties (e.g., geometry and depth) and magnitude of the imposed overpressure (ΔP). For example, a model incorporating a high reservoir temperature with a low geothermal gradient (e.g., 1100°C and 20 K km −1 ), which produces increased long-term displacements (Figure 9b), may underestimate the magnitude of overpressure or increase the depth of the best-fit source relative to a model that couples a low reservoir temperature with a high geothermal gradient (e.g., 700°C and 50 K km −1 ). Whilst the use of a temperature-dependent inelastic rheology likely provides an improved understanding of the source processes underlying an episode of unrest, relative to elastic models, this analysis highlights the potential uncertainty that can arise if the imposed thermal regime is poorly constrained. In this study, we consider the application of constant and ramped loading, however, the impact of uncertainties in the thermal regime will also be observed with alternative temporal loading functions, such as rate-decreasing, trapezoidal or multiramped loading functions, due to the Boltzmann superposition principle.

Comparisons to Viscoelastic Shell Models
A common implementation of Maxwell viscoelasticity in ground deformation studies is as a shell surrounding the deformation source, embedded within an elastic half-space (e.g., Delgado et al., 2018;Dragoni & Magnanensi, 1989;Newman et al., 2001;Novoa et al., 2019), as illustrated in Figure 10a. This approach assumes that viscoelastic behavior is restricted to the thermal aureole surrounding a magmatic system whilst the surrounding crust behaves in an elastic manner, allowing deviatoric stresses to be maintained outside the shell and modulating the ductile response. Upon application of a constant overpressure, this configuration produces an instantaneous deformation response, according to the source radius, 1 R , followed by rate-decreasing displacement bounded by the elastic response of an effective source with a radius of the viscoelastic shell, 2 R (Bonafede & Ferrari, 2009;Dragoni & Magnanensi, 1989;Segall, 2010). The characteristic timescale of the viscous deformation is primarily controlled by the shell viscosity, , but is also dependent on the geometry of both the source and shell (Figure 10a). Newman et al., (2001) suggest that the timescale of rate-decreasing displacement is a useful parameter to describe the deformation episode, which could be used to estimate the properties of the viscoelastic shell based on goodness-of-fit (e.g., Delgado et al., 2018; HEAD ET AL.
The Maxwell viscoelastic models of our study can be considered as an endmember of the viscoelastic shell modeling approach, whereby the shell radius, 2 R , is no longer finite. This presents an opportunity for critical analysis of using thermally informed viscoelastic shells against a thermo-viscoelastic model. Using the viscosity distribution for a reservoir temperature of 900°C and a geothermal gradient of 30 K km −1 , the geometry for a series of discrete, nonuniform viscoelastic domains is determined from the viscosity contours shown in Figure 10b. The minimum thickness of each shell occurs at the top of the modeled reservoir, and is approximately 340, 620, 900, and 1,160 m thick for the 10 16 , 10 17 , 10 18 , and 10 19 Pa s contours, respectively. Based on the geothermal gradient used, the 10 19 Pa s contour extends to the lateral edge of the model space and therefore captures the viscosity structure to the base of the model. These shells are of heterogeneous viscosity, as determined by the temperature distribution, and are embedded within an elastic domain. Comparing the resultant deformation timeseries, as shown in Figure 10c, highlights the potential disparity in modeled displacements when using discrete viscoelastic shells. Restricting the geometry of the viscoelastic domain decreases the magnitude of the resultant uplift and reduces the deformation timescale (e.g., Bonafede & Ferrari, 2009). Using a maximum viscosity of 10 19 Pa s provides a close match to the magnitude of displacement produced by the thermo-viscoelastic model, but this is largely due to the ∼annual-scale to ∼decadal-scale timescales captured, and the modeling period considered. Due to the difference in displacements produced by the thermo-viscoelastic model and those with a viscoelastic shell (Figure 10c), it is clear that modeling a deformation episode with a restricted viscoelastic geometry, for example, the <10 16 Pa s shell model, will require increased overpressures, or a shallower deformation source, compared to the thermo-viscoelastic model. Whilst there are rheological arguments against the use of the Maxwell viscoelasticity (characteristics of fluid behavior, including limitless creep and stress relaxation; e.g., Banks et al., 2011;Marques & Creus, 2012), this analysis presents a further complication for using a viscoelastic shell approach when modeling volcanic ground deformation; raising the question as to whether the best-fit underlying source and inferred temporal loading are accurate.
HEAD ET AL.

Limitations of the Modeling Approach
The model structure used in this study assumes a discrete magmatic reservoir emplaced at a mid-crustal depth. As the reservoir is represented by a cavity, our models do not account for internal magma dynamics, both as driving forces and second-order effects. Integrating magma dynamics within deformation models would provide valuable insight into the possible mechanisms that can produce the pure pressure-change and volume-change deformation modes explored here, or whether a hybrid of the two modes is likely to be more prevalent. Furthermore, we use an idealized reservoir geometry, whereas the aspect ratio of deformation sources varies widely across studies of volcanic ground deformation. Whilst this does not raise an issue when comparing the results for a single geometry, as is the case in this investigation, the aspect ratio of the reservoir is a first-order control on the partitioning of the deformation field. As a result, oblate geometries, which preferentially partition deformation vertically (if they are not dipping), may be more greatly affected by the partitioning that arises from viscosity gradients beneath the reservoir. Furthermore, whilst not explicitly modeled, the presence of a vertically extensive transcrustal magmatic system (e.g., Cashman et al., 2017) beneath the reservoir would likely perturb the regional geothermal gradient to greater depths, due to the thermal priming and maturation of the crustal section (Annen, 2011;Karakas et al., 2017). A future consideration could be to incorporate lateral temperature gradients, such that perturbations to the background geotherm are more localized, to investigate how the presence of a transcrustal magmatic system would impact the partitioning of deformation and the resultant surface displacements. Furthermore, a mush-dominated transcrustal magmatic system may accommodate an increased component of deformation beneath the reservoir, which likely necessitate localized regions with alternative behaviors such as poroelasticity (e.g., Liao et al., 2018Liao et al., , 2021. The results of this study outline the importance of the thermal constraints used, and therefore the resultant viscosity distribution, when implementing thermo-viscoelasticity in ground deformation models. Ultimately, this highlights the assumed relationship between the steady-state temperature field and the viscosity; the Arrhenius formulation (Equation 2). This states that the viscosity increases exponentially, inversely proportional to the temperature field. However, the resultant viscosity is also controlled by the Dorn parameter, a pre-exponential factor, and the activation energy of the surrounding host rock. An investigation into the parameter space of the Arrhenius formulation by Morales Rivera et al. (2019) demonstrates that the Dorn parameter and activation energy can be altered to account for different crustal compositions (Kirby & Kronenberg, 1987;Ranalli & Murphy, 1987), with the capacity to alter the resultant viscosity distribution by several orders of magnitude. However, caution should be exercised when choosing parameter values, as the availability of data relating to different rock compositions is limited and the scaling of parameters from laboratory studies to crustal-scale applications is uncertain. In this study, we use felsic parameters that are employed in similar investigations (Del Negro et al., 2009;Gregg et al., 2012;Hickey et al., 2016). As such, in order to improve the representation of crustal viscosity in deformation models, it is important to understand how values of the Dorn parameter and activation energy may be expected to change in a volcanic setting (i.e., with composition, depth/confining pressure, temperature, and alteration), and to combine complementary data sets to determine the local thermal structure.
Both the viscoelastic rheologies explored in this study are highly idealized, whereby simple relationships between the spring and dashpot components (Figures 1b and 1c) produce nonelastic rheological effects. While we make the recommendation for the use of the SLS rheology over Maxwell, it is worth noting that it still unlikely provides a complete representation for the rheology of the entire crustal column, whereby the ductility of the lower crust accommodates finite strain over long (geological) timescales. However, it is more robust than the Maxwell model, and offers simplicity (in terms of undefined parameters) over more complex power-law rheologies, where the effective viscosity depends on the imposed strain-rate (e.g., Pearse & Fialko, 2010;Ranalli, 1995). Furthermore, a viscoelastic representation of the crust does not account for components of plasticity, where, in reality, a finite yield threshold can maintain far-field deviatoric stresses. Understanding the rheology of crustal rocks, through laboratory and field-based experiments, as well as the scope and parameter-space of rheological models, in conceptual studies such as this, are critical for identifying the mechanisms driving ground deformation, as well as a broad spectrum of Earth science problems.

Conclusions
Our study indicates that the choice of the modeling approach (i.e., viscoelastic rheology, deformation mode, and chosen thermal regime) has significant implications for interpreting geodetic volcano monitoring data. The observed deformation is a function of the magmatic processes and the rheological response of the surrounding crustal material, both of which can vary through time. By considering a range of thermal regimes, we demonstrate that thermomechanical heterogeneity results in markedly different surface displacements, despite an identical underlying mode of deformation. Whilst many studies of volcanic deformation consider the heterogeneity of elastic moduli within the crust, we highlight that the thermal regime, and therefore viscosity structure, plays a primary role in affecting the predicted ground deformation for a given viscoelastic rheology. Usually, viscous effects are thought to be limited to long-term deformation responses, and negligible over the short-term, however the viscosity structures considered in this study demonstrate that thermo-viscoelasticity is capable of capturing a range of deformation timescales in response to an imposed overpressure (ΔP) load. High rates of deformation over short timescales can be produced, provided that low viscosities relating to high-temperature crustal material in the vicinity of the modeled reservoir, that generate the short-timescale signals, are included. Long-term deformation signals, such as steady uplift, may not represent ongoing activity within the magmatic system, but instead might reflect the response of cooler crustal materials, of greater viscosity and longer timescales, closer to the ground surface. This may be compounded by variations in pore pressure and the thermal expansion of crustal rocks (e.g., Coco et al., 2016).
Comparisons of the Maxwell and SLS rheological models demonstrate the disparity between the deformation responses of a viscoelastic fluid and a viscoelastic solid in the presence of viscosity gradients, with the fluid behavior of the Maxwell rheology partitioning elevated deformation into viscous regions and producing a long-term subsidence signal. Whilst Maxwell viscoelasticity is often used in deformation modeling when investigating inelastic effects due to its simplicity, its variable behavior within models that are currently best thought to represent reality (i.e., heterogeneous crustal viscosity), arising from intrinsic fluid behavior, raises questions about the applicability of the rheology to models of volcanic ground deformation. Given that the SLS rheology exhibits a consistent, finite deformation response, we recommend that its use should be favored. Incorporating viscoelasticity generally acts to reduce the overpressure requirements suggested by elastic modeling, due to the increased magnitude of the deformation response. Here, we identify that uncertainties in the imposed thermal constraints (i.e., reservoir temperature and geothermal gradient) may appreciably impact the inferences of the best-fit source parameters (e.g., size and depth) and overpressure loading history from thermo-viscoelastic models, as the resultant viscosity distribution affects deformation partitioning.

Data Availability Statement
Data were not used nor created for this research.