The Influence of Viscoelastic Crustal Rheologies on Volcanic Ground Deformation: Insights From Models of Pressure and Volume Change

Inelastic rheological behavior, such as viscoelasticity, is increasingly utilized in the modeling of volcanic ground deformation, as elevated thermal regimes induced by magmatic systems may necessitate the use of a mechanical model containing a component of time‐dependent viscous behavior. For the modeling of a given amplitude and footprint of ground deformation, incorporating a viscoelastic regime has been shown to reduce the magma reservoir overpressure requirements suggested by elastic models. This phenomenon, however, is restricted to pressure‐based analyses and the associated creep behavior. Viscoelastic materials exhibit additional constitutive time‐dependent behaviors, determined by the stress and strain states, that are yet to be analyzed in the context of volcanic ground deformation. By utilizing a mechanically homogeneous model space and distinct reservoir evolutions, we provide a comparison of three viscoelastic rheological models, including the commonly implemented Maxwell and Standard Linear Solid configurations, and their time‐dependent behaviors from a fundamental perspective. We also investigate the differences between deformation time series resulting from a pressurization or volume change, two contrasting approaches that are assumed to be equivalent through elastic modeling. Our results illustrate that the perceived influence of viscoelasticity is dependent on the mode of deformation, with stress‐based pressurization models imparting enhanced deformation relative to the elastic models, thus reducing pressure requirements. Strain‐based volumetric models, however, exhibit reduced levels of deformation and may produce episodes of apparent ground subsidence induced by source inflation or vice versa, due to the relaxation of crustal stresses, dependent on whether the reservoir is modeled to be expanding or contracting, respectively.


Introduction
Distinguishing the underlying processes that drive episodes of volcanic unrest is of great importance for understanding the behavior of a subvolcanic system and elucidating its evolution. The measurement of ground deformation, among other observables, in active volcanic regions can provide insights into the mechanisms driving unrest, such as the migration and accumulation of magma (e.g., Bato et al., 2018), the cooling and crystallization of magma (e.g., Dzurisin et al., 1990) and exsolution of volatiles ©2019. The Authors. This is an open access article under the terms of the Creative Commons Attribution License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.

10.1029/2019JB017832
Key Points: • The influence of viscoelasticity depends on whether it is utilized in pressurization (stress-based) or volume change (strain-based) models • Time-dependent creep and recovery behaviors of the popular Maxwell viscoelastic configuration can result in large irreversible strains • More realistic pressure requirements can be derived from elastic models using scaling factors independent of source and material properties (e.g., Caricchi et al., 2014), or changes to hydrothermal systems (e.g., Fournier & Chardot, 2012). Identifying the processes that are responsible for the observed deformation, in turn, enhance the knowledge of the underlying magmatic plumbing system and the hazards posed by unrest episodes (Sparks, 2003).
Traditionally, the observed deformation field is modeled analytically, as a point source (Mogi, 1958) within an isotropic and homogeneous elastic half-space, providing a very simple interpretation of the sources responsible for the observed deformation. Over time, these models have evolved to account for a variety of finite source geometries (Fialko et al., 2001;McTigue, 1987;Okada, 1985;Yang et al., 1988), as well as accommodating additional complexities, including spatially variable components such as topography (e.g., Cayol & Cornet, 1998), structural discontinuities (e.g., De Natale et al., 1997), and medium heterogeneity (e.g., Trasatti et al., 2003), with a shift toward numerical modeling. While elastic models are widely utilized and are capable of reproducing uplift and subsidence patterns observed at volcanic centers, they often rely on pressure changes with unrealistic amplitudes (e.g., Del Negro et al., 2009;Masterlark et al., 2010;Newman et al., 2006). The use of an elastic rheology presents a simple foundation for deformation modeling; however, this approximation is generally only appropriate for the deformation of crustal materials at temperatures cooler than the brittle-ductile transition, which ranges from 300 to 600°C dependent on the strain rate and crustal composition (Del Negro et al., 2009;Newman et al., 2001;Ranalli, 1995), which occur over short timescales and result in small displacements. Hence, additional rheological effects can be incorporated when modeling observed deformation patterns, especially in volcanic regions where perturbations to the crustal geotherm, induced by long-lived magmatic systems (e.g., Annen, 2011;Gelman et al., 2013;Karakas et al., 2017), are expected to characterize the behavior of the middle crust and upper crust (e.g., de Silva & Gregg, 2014). Elevated thermal regimes allow crustal materials to behave in a nonelastic manner, where a rheological model containing a time-dependent (viscous) component of behavior is more likely to represent the observed deformation field (Jellinek & DePaolo, 2003;Newman et al., 2001), despite the apparent goodness-of-fit provided by elastic models. Early models that incorporated viscoelastic rheologies were limited to the allocation of a single viscosity across the model space (Bonafede et al., 1986;Segall, 2010) or consisted of a viscoelastic shell encompassing the source of deformation (Bonafede & Ferrari, 2009;Dragoni & Magnanensi, 1989). A fundamental characteristic of viscoelasticity is the ability to more accurately represent deformation time series, due to intrinsic time-dependent behaviors. Elastic models coupled with other phenomena, such as flow within magmatic plumbing systems (e.g., Le Mével et al., 2016;Lengliné et al., 2008), are also capable of producing time-dependent deformation signals proximal to a volcanic centre. Due to the strong influence of temperature on viscosity, thermal profiles have been incorporated into viscoelastic models to account for heterogeneous rheological properties, which are known to affect the partitioning of stress and strain in volcanic settings (Del Negro et al., 2009;Gottsmann & Odbert, 2014;Gregg et al., 2012;Hickey et al., 2015;Hickey et al., 2016).
With a variety of different deformation models available both analytically and numerically, it is important to ensure that the choice of model parameters provides an accurate representation of the modeled region. If the rheology utilized for a model is unrealistic, despite the construction of numerous models that will satisfy a goodness-of-fit criterion to a given deformation time series, the inferences of the underlying deformation processes are unlikely to be correct. As such, it is important to have a fundamental understanding of how different viscoelastic rheologies behave under different states of stress and strain within simple homogeneous model spaces, and how their time-dependent behaviors are affected, prior to the construction of models with greater degrees of complexity. Here, we investigate the influence of different viscoelastic crustal rheologies on ground deformation resulting from a subsurface deformation source, representing a magma reservoir, which is allocated either a pressurization or volume change deformation mode. Analytical modeling techniques often resolve volumetric changes when identifying episodes of volcanic unrest (e.g., Parker et al., 2014), whereas numerical methods often use changes in overpressure to reproduce deformation time series (e.g., Hickey et al., 2016). The difference between these source deformation modes, stress-based (ΔP) versus strain-based (Δα), in a viscoelastic regime is often overlooked, as traditional elastic models generally consider them to be equivalent. However, when incorporating a rheology that has behaviors dependent on the stress and strain states, these deformation modes impart significant differences in the modeled spatial and temporal deformation patterns and hence the inferences of a magmatic unrest episode.

Linear Viscoelastic Model
Hooke's law is the fundamental principle underlying elasticity, stating that the strain, ε, in a solid is proportional to the applied stress, σ. The uniaxial (1-D) relation takes the form: where the coefficient of proportionality is the elastic modulus, E. Deformation is defined as being elastic when the induced strains occur instantaneously and can be fully recovered, resulting in the material returning to its original form as the strains are removed, and the deformation itself is independent of time (Christensen, 1982;Dzurisin, 2007). However, materials cannot be considered to behave elastically if the relationship between stress and strain is variable with time. In this scenario, viscous effects must be considered, whereby the strain rate, _ ε, is directly proportional to the applied stress, σ: where the coefficient of proportionality is the viscosity, η. In reality, very few materials exhibit perfect elastic behaviors beyond small deformations and very short deformation timescales (Christensen, 1982;Ranalli, 1995). As a result, many materials demonstrate a combination of both instantaneous elastic and timedependent viscous behaviors as the magnitude and timescale of deformation varies, and so are considered to be viscoelastic (Christensen, 1982;Crawford, 1998;Ranalli, 1995). Viscoelastic rheologies are increasingly utilized in the modeling of crustal deformation, which is a particularly important consideration in volcanic settings. Elevated thermal regimes surrounding magmatic bodies are believed to invalidate the use of an elastic approximation, due to increasing temperatures raising the brittle-ductile transition to shallower levels and increasing the prominence of a viscous behavioral component (Dragoni & Magnanensi, 1989;Newman et al., 2001). A common assumption when implementing viscoelasticity is that the viscous component of deformation is incompressible, which allows for volumetric strains to be considered as purely elastic and the viscoelastic deformation to be represented in terms of the deviatoric components (Segall, 2010). As such, the bulk modulus, K, of the material behaves elastically, whereas the shear modulus, G, behaves in a viscoelastic manner (Del Negro et al., 2009;Folch et al., 2000). Several simple viscoelastic configurations can be conceptually represented by the linear combination of elastic springs, with a spring constant of shear modulus G, and viscous dashpots, with a viscosity coefficient of η, which govern the contributions of the elastic and viscous components to the response of the rheological models to changes in stress and strain ( Figure 1).
Viscoelastic materials exhibit three fundamental time-dependent behaviors, creep, relaxation, and recovery, which are each dependent on the states of stress and strain that have been applied to the material ( Figure 1). Creep behavior occurs when a material is subjected to a constant stress ( dσ dt ¼ 0 ) and describes the associated increase in strain (Christensen, 1982;Crawford, 1998). If instead a material is subjected to a constant strain ( dε dt ¼ 0), then the material may exhibit a dissipation of stress, a behavior known as relaxation (Christensen, 1982;Crawford, 1998). Upon the removal of an external stress, a material undergoes recovery, which describes the removal of strain within the material, and is analogous to the creep behavior (Christensen, 1982;Crawford, 1998). In this investigation, we consider the behaviors of linear viscoelastic configurations, in which stress depends linearly on the strain and strain rate, and we include the Maxwell, Kelvin-Voigt (KV), and Standard Linear Solid (SLS) rheological models in our analysis.

Maxwell
The simplest viscoelastic configuration is the Maxwell model, formed of a spring and dashpot in series as shown in Figure 1a. This arrangement results in stress, σ M , being applied equally across both elements, while the total strain, ε M , is the sum of the contributions from each component (Christensen, 1982;Crawford, 1998), given below: where _ ε M and _ σ M denote the time derivatives of the strain and stress, respectively. 2.1.1. Creep Following the application of a constant stress, σ A , the strain as a function of time is described by the following (Christensen, 1982;Crawford, 1998): This equation illustrates that the creep behavior is formed of two components, an instantaneous, elastic strain, σA E , and a viscous strain, σA η t, which is a linear function of time.

Recovery
If the strain is held constant, and so _ ε M ¼ 0, the stress as a function of time is given by (Christensen, 1982;Crawford, 1998) where σ 0 is the stress at the onset of constant strain. This equation indicates that the stress relaxation behavior features an exponential decay with a time constant of η E , referred to as the relaxation, or Maxwell, time. If the rheology is allowed to relax for sufficiently long period of time, it will approach a state of zero stress.

Relaxation
If a stress, σ′, is removed from the system, the level of strain is given by  Crawford (1998). The time-dependent behavior schematics illustrate creep behavior under a constant stress ( dσ dt ¼ 0), the relaxation response to a constant strain ( dε dt ¼ 0), and the recovery of strain following the removal of stress. The elastic strain (denoted ε E ) imparted by a constant stress is represented by a dashed horizontal line within the strain plots.
where ε 0 is the strain prior to stress removal. There is an instantaneous recovery of the elastic strain, ε 0 − σ 0 E , proportional to the stress removed. Following this, the strain rate is then zero and there is no further viscous strain recovery.

Summary
These behaviors are summarized within Table 1 and schematically in Figure 1a. It can be seen that the Maxwell viscoelastic model provides a first-order approximation of the viscous stress relaxation when subjected to a constant strain; however, it does not accurately depict the expected creep or recovery behaviors of crustal materials. This is a result of the strain, and hence deformation, increasing linearly without bound under the application of a constant stress, and a permanently deformed state following the removal of the initial stress.

Kelvin-Voigt
The Kelvin-Voigt (KV) rheological model, as shown in Figure 1b, is formed of a spring and dashpot in parallel. This configuration, in contrast to the Maxwell model, subjects each element to the same strain, ε KV , whereas the total stress within the system, σ KV , is the sum of the contributions from each component, given below: Using the above criteria, and the constitutive relations (1) and (2), the governing equation for the KV model (Christensen, 1982;Crawford, 1998) is given by where _ ε KV is the strain rate.

Creep
The strain as a function of time, following the application of a constant stress, σ A , is described by the equation (Christensen, 1982;Crawford, 1998): Accordingly, the creep behavior of this model results in the exponential strain increase from zero to the asymptote σA E , the value of the elastic strain response, with a time constant of E η . As there is no strain Relaxation (constant strain) Stress decays exponentially No relaxation, constant stress Stress decays exponentially Recovery (removal of stress) Instant recovery of elastic strain, no viscous recovery Exponential recovery of strain Instant then exponential recovery of strain Note. σ A denotes the applied stress, σ′ indicates the stress removed, while σ 0 and ε 0 represent the prior states of stress and strain, respectively. The relative weights of the Standard Linear Solid arms are given by μ 0 and μ 1 , corresponding to the spring and Maxwell arms respectively.
response with the application of stress at time t = 0, it is apparent that the KV rheology lacks the instantaneous elastic strain observed in the Maxwell rheology.

Relaxation
If a constant strain is applied to the model, then the governing equation decomposes to a form that is independent of time, and so the stress remains constant. The model is then supported by the spring element, producing the response of an elastic material that does not exhibit any stress relaxation behavior.

Recovery
If the stress is removed from the system, and so σ = 0, the strain as a function of time is described by the following (Christensen, 1982;Crawford, 1998): where ε 0 is the strain at the time of stress removal. This equation represents an exponential recovery of the strain, with a time constant of η E , which is a reversal of the predicted creep for this rheological model.

Summary
These behaviors are summarized within Table 1 and schematically in Figure 1b. A strength of the KV model is that the creep behavior is asymptotic in nature, which could be considered a weakness of the Maxwell model, although it is bounded by the elastic response of the spring. Despite this, the KV model is likely unsuitable for the representation of many materials, including crustal rocks, as it cannot produce an instantaneous elastic strain response when subjected to a load. Further to this, this model is unable to account for the relaxation of stress when strained, and an instantaneous strain application results in an infinite stress response.

Standard Linear Solid
The Standard Linear Solid (SLS) model, as seen in Figure 1c, is produced by combining a spring and a Maxwell arm in a parallel configuration. As a result, the total stress of the system, σ SLS , is given by the sum of the stresses across both the Maxwell and spring arms, whereas the strain is equal in both arms and so the total strain, ε SLS , is then given by the sum of the strains within the Maxwell arm. This is illustrated below: Utilizing these criteria and their time derivatives, the governing equation for the SLS rheological model takes the form (Christensen, 1982;Crawford, 1998) where μ 0 and μ 1 are the fractional shear moduli in the elastic and Maxwell arms, respectively. As a result, the end-member solutions are either elastic or Maxwell in nature; however, the specific relative contribution for a given material is an unknown factor in the modeling process. In this work we consider both arms to be equally weighted, where μ 0 = μ 1 = 0.5.

Creep
With the application of a constant stress, σ A , the strain as a function of time is given by (Christensen, 1982;Crawford, 1998) The creep behavior is composed of the instantaneous elastic strain, σA E , followed by an exponential strain increase with a time constant of η μ 0 μ 1 E . If the applied stress is kept constant for a sufficient period of time, the strain will approach a value of σA μ 1 E , which is controlled by the relative contribution of the Maxwell arm.

Relaxation
If the strain, however, is held constant, _ ε SLS ¼ 0, the resultant stress at any given time is described by (Christensen, 1982;Crawford, 1998) where σ 0 is the level of stress prior to the constant strain. This states that the stress within the system will decay exponentially with a relaxation time of η μ 1 E , which is controlled by the relative contribution of the Maxwell arm. If the system is permitted to undergo relaxation for a sufficient amount of time, the stress will asymptote to a value, μ 0 σ 0 , that is controlled by the relative contribution of the elastic arm.

Recovery
If a stress, σ′, is removed from the system, the time-dependent strain can then be determined by (Christensen, 1982;Crawford, 1998) where ε 0 is the strain prior to stress removal. This predicts an instantaneous recovery of the elastic strain ε 0 ¼ σ 0 E , proportional to the amount of stress removed from the system. This is then followed by an exponential decay of the 'viscous' strain, controlled by a time constant of η μ 0 μ 1 E , with the strain returning to zero if allowed to recover over a sufficient timescale.

Relaxation
These behaviors are summarized within Table 1 and schematically in Figure 1c. The SLS rheological configuration is the simplest combination of components to provide a reasonable representation of each time-dependent viscoelastic behavior, by including creep behavior that contains both an instantaneous (elastic) and an asymptotic timedependent (viscous) component, an asymptotic stress relaxation response, and the full recovery of strain in the absence of stress.

Rheological Considerations
With the significant differences in the behaviors for each of the viscoelastic configurations described above, it is clear that the inferences drawn from the temporal evolution of a magmatic system undergoing unrest are highly dependent on the modeled rheology. The prominent characteristics for each of the viscoelastic configurations and the associated time-dependent behaviors are summarized within Table 1, while full derivations and symbols can be found in the Supporting Information.

Setup and Geometry
We use COMSOL Multiphysics® (v5.3a) to construct and solve finite element (FE) forward models of ground deformation with the structural mechanics module. We utilize a 2D-axisymmetric model geometry containing homogeneous half-space characteristics, which are representative of those found in volcanic regions (Gudmundsson, 2011), and simple source parameters (Table 2). We do not explicitly account for the effect of gravitational loading within our models in order to reduce computational cost, as this only affects the distribution of the resultant stress, and hence has implications for the locations of reservoir failure (e.g., Grosfils, 2007), but does not the influence the overall pattern of deformation. By employing a modeling approach with few complexities, we ensure  that the study focuses foremost on the differences between the rheological models from a fundamental perspective.
The model space is adapted from the benchmarked approach outlined in Hickey and Gottsmann (2014), represented schematically in Figure 2, consisting of vertical and horizontal dimensions of 52 km. The deforming magma reservoir is represented by a finite spherical cavity of radius (α) 1,000 m, which is centered at a depth (d) of 5,000 m below the free surface. In parallel model setups, the boundary of the cavity is allocated one of two different deformation modes, either an overpressure (ΔP) or a prescribed expansion (Δα). In order to ensure agreement between these deformation modes, the degree of expansion is derived from elastic models, with a source overpressure of 10 MPa, by evaluating the resultant volume change along the reservoir boundary within COMSOL Multiphysics®. The volume change (ΔV) owing to this pressurization is calculated to be~3.94 × 10 6 m 3 , providing an equivalent, radially uniform, expansion of~31.36 cm that is applied to the source boundary. In this investigation, we consider three different half-space viscosities (10 17 , 10 18 , and 10 19 Pa s) and evaluate the vertical deformation time series directly above the center of the source, at z = r = 0, which are normalized to the response of the corresponding elastic models. Time series that have not been normalized can be found in Supporting information Figures S2-S5. Both the elastic and viscoelastic models are executed with time-dependent deformation modes, over a time period of 10 years with a temporal resolution of 0.01 year.

Temporal Reservoir Evolutions
In order to provide an extensive comparison of each of the viscoelastic behaviors, we consider four distinct time evolutions for the deformation modes of the modeled magma reservoir, the magnitude of which evolves over a time period of 10 years ( Figure 3). This chosen time period is of sufficient length to allow the time-dependent behaviors of each viscoelastic configuration to be conveyed within the range of modeled parameters and is relevant to volcano monitoring timeframes (Phillipson et al., 2013). The deformation modes implemented in this investigation are static and so neglect any dynamical pressure-volume relationships, such as the reduction in overpressure following inflation (e.g., Gregg et al., 2013). The four functions used in our analysis are illustrated in Figure 3 and include constant, ramped, rectangular, and trapezoidal inputs, which are discussed in their respective results sections. Traditionally with elastic modeling, the inferred temporal reservoir evolution, and hence overpressure and expansion requirements, directly reflects the profile of a deformation time series.

Results
Here we compare the resultant vertical deformation time series for 24 different models, formed by a combination of the different reservoir evolutions (4), viscoelastic rheologies (3), and the modes of deformation (2). Time-independent elastic models are evaluated for each reservoir evolution and deformation mode, and each viscoelastic model is tested over three half-space viscosity values (10 17 , 10 18 , and 10 19 Pa s) to further demonstrate the viscoelastic behaviors. The fundamental differences between the stress-based (ΔP) and strain-based (Δα) deformation modes are driven by the creep and relaxation responses to changes in stress and strain conditions, respectively. Creep behavior describes the time-dependent increase in strain, in response to a constant stress, and so elevates the level of observed deformation, whereas the relaxation response describes the time-dependent dissipation of stress in the presence of a constant strain, and so reduces the observed deformation over time. The analysis of results below is broken down by first considering the reservoir evolution, followed by the deformation mode, and then the implemented rheology.

Constant Forcing
The constant reservoir condition is displayed in Figure 3a and may 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. The resultant model outputs are presented in Figure 4, with pressurization and expansion model results occupying the upper and lower rows respectively. With this time function, the pressurization deformation mode exhibits a constant stress, and so we observe creep behavior, whereas the expansion deformation mode (equivalent to a volume change) emplaces a constant strain, and so relaxation behavior is recognized instead.

Pressurization
Upon first inspection of the pressurization models, both the Maxwell and SLS models result in amplified vertical deformation relative to the elastic model, the magnitude of which is dependent on the viscosity of the model space, whereas the KV models exhibit a convergence to the elastic solution. With the Maxwell models we observe exaggerated deformation up to 18 times greater than the elastic model over a time period of 10 years, for the lowest modeled viscosity, which is a result of the unbounded linear creep behavior. In contrast, the SLS models demonstrate deformation up to a factor of~1.7 times greater than the elastic model, which constitutes an asymptote for the 10 17 Pa s viscosity model. The KV models are distinct from the Maxwell and SLS models due to higher viscosity values producing a greater deviation from the elastic model. This is attributed to a combination of the rheology lacking an elastic strain response and the eventual convergence to the elastic solution, a consequence of the reduced rate of deformation with increasing viscosity.

Expansion
With the expansion deformation mode, the reservoir boundary undergoes an instant expansion at t = 0, and so the system is immediately strained. We observe for both the Maxwell and SLS models that they produce the same magnitude of deformation as the elastic model upon expansion. Following this, these models exhibit subsidence behavior due to the dissipation of stress, a result of the relaxation behavior in response to a constant strain. It can be seen that the deformation time series for both the Maxwell and SLS rheologies converge to a permanently uplifted state of~0.66 and~0.85 times that of the elastic deformation, respectively, the rate of which is dependent on the viscosity of the model space. In contrast, deformation for  (Figure 3a) reservoir condition, with the time series normalized to the elastic solutions. By comparing the top row with the bottom row, the time series illustrate fundamental differences between the stress-based (ΔP) and strainbased (Δα) deformation modes, and the associated creep and relaxation responses, for each viscoelastic rheology. Differences between the responses to the same deformation mode (i.e., ΔP or Δα, across a row), for the Maxwell, Kelvin-Voigt, and Standard Linear Solid models, highlight the importance of carefully considering a chosen rheology. The elastic vertical deformation is 3.76 and 3.11 cm for the pressurization and expansion deformation modes, respectively.

10.1029/2019JB017832
Journal of Geophysical Research: Solid Earth the KV models originates from zero rather than that of the elastic model. This is thought to arise from a numerical simplification within COMSOL Multiphysics®, as the constitutive equations for the KV model state that an instantaneous change in strain results in an infinite stress (Marques & Creus, 2012). Following the strain change at t = 0, we observe that the KV deformation time series again converge to the elastic solution, as seen in the pressurization models. As this rheology does not allow for stress relaxation, an elastic solution is maintained.

Ramped Forcing
The ramped reservoir evolution, shown in Figure 3b, is formed of a single linear segment and may represent the sustained, constant-rate injection of magma or exsolution of volatiles. The corresponding model results are displayed in Figure 5, with pressurization and expansion model results occupying the upper and lower rows. With this time function, the pressurization deformation mode exhibits a linearly increasing stress, whereas the expansion deformation mode emplaces a constant strain rate. The resultant time series are controlled by the Boltzmann superposition principle, which states that each increment of load makes an independent and additive contribution to the total deformation (Vincent, 1982), and so we observe compounded creep and relaxation behaviors with each time step. In a similar fashion, compounded behaviors would also be observed for variations on the linear ramp input, such as exponential or logarithmic reservoir evolutions.

Pressurization
In accordance with the Boltzmann superposition principle, we observe Maxwell and SLS deformation time series that are amplified relative to elastic solution, a result of the compounded creep behavior. With the increase in stress that occurs with each time step, there is an independent strain contribution to the total deformation that is controlled by the creep behavior of the viscoelastic model. Specifically, for the Maxwell model, the linearly unbounded creep response results in exponentially increasing levels of deformation. The asymptotic nature of the SLS model (Figure 4), however, provides a limit to the rate of deformation that is produced, as the rate of deformation becomes linear once the cumulative strains reach the asymptote. This occurs at around t = 8 years for the 10 17 Pa s model. As the creep behavior of the KV The results demonstrate fundamental differences between the stress-based (ΔP) and strain-based (Δα) deformation modes for each viscoelastic rheology, and the associated compound creep and relaxation responses, by comparing the top row with the bottom row. The differences between the responses of the Maxwell, Kelvin-Voigt, and Standard Linear Solid models, for the same deformation mode (i.e., ΔP or Δα, across a row), highlight the importance of carefully considering a chosen rheology. The maximum elastic vertical deformation is 3.76 and 3.11 cm for the pressurization and expansion deformation modes, respectively. rheological model lacks an instantaneous elastic strain component, reduced levels of deformation are observed relative to the elastic solution.

Expansion
In contrast to the pressurization models, the deformation resulting from the Maxwell and SLS expansion models do not exceed the elastic solution. With each time step, the same amount of strain is applied to the crust surrounding the expanding reservoir, resulting in uniform increases in deformation. The deviation we observe between the Maxwell and SLS models and the elastic solution, however, are due to the viscous stress relaxation that occurs between the strain-additive time steps. This enables the strain-induced stress response to partially dissipate before it can contribute to the observed inflation. As the rate of relaxation is viscosity-dependent, more relaxation occurs between each time step for the lower viscosity models, and therefore, the deformation time series increases at a slower rate, resulting in a greater deviation from the elastic solution. Due to the asymptotic nature of the stress relaxation response (Figure 4), the long-term deformation rate for both the Maxwell and SLS models becomes linear. The KV models appear to follow this same process; however, greater deformation is exhibited for higher-viscosity values, in contrast to the Maxwell and SLS models. This is due to an exaggerated stress response with each strain increase, which is believed to be a numerical simplification within COMSOL Multiphysics®, as the stress response should be infinite for instant strain changes (Marques & Creus, 2012). Between each time step, prior to the addition of an additional strain, the existing strain undergoes a decay toward the elastic solution, with a viscositydependent rate. This results in greater viscosities decaying more slowly, producing enhanced deformation with the addition of another strain. This effect is readily observed in section 4.3.2.

Rectangular Forcing
The rectangular reservoir evolution (Figure 3c), formed of a unit step (Heaviside) function and its subsequent reversal, describes the instant accumulation of magma or the exsolution of volatiles at a time of 3 years, followed by an instant dissipation at 7 years, which may be in the form of an intrusion or eruption. Between these two stages there is a 4-year period in which the reservoir boundary undergoes no change, which may  (Figure 3c) reservoir evolution, with the time series normalized to the elastic solutions and the elastic strain denoted by ε E . Fundamental differences between the stress-based (ΔP) and strain-based (Δα) deformation modes are observed for each viscoelastic rheology, by comparing the top row with the bottom row. The pressurization models exhibit creep and recovery behaviors, in response to the constant reservoir condition and the removal of forcing, respectively, whereas the expansion models undergo episodes of relaxation. Further differences for the same deformation mode (i.e., ΔP or Δα, across a row), between the Maxwell, Kelvin-Voigt, and Standard Linear Solid models, highlight the importance of carefully considering a chosen rheology. The maximum elastic vertical deformation is 3.76 and 3.11 cm for the pressurization and expansion deformation modes, respectively. reflect a state of equilibrium. The model results are displayed in Figure 6, with pressurization and expansion model results occupying the upper and lower rows respectively. Fundamental differences between the pressurization and expansion models are observed with the constant reservoir condition imposed between t = 3 years and t = 7 years, due to the respective creep and relaxation responses. Another key difference is observed at t = 7 years, with pressurization models exhibiting strain recovery behavior, following the removal of stress, whereas the expansion models undergo stress relaxation due to an effective negative strain.

Pressurization
We observe each of the rheological models exhibiting the same creep behavior as seen in the constant overpressure models (Figure 4) at the onset of pressurization at t = 3 years. While the timescale of the constant overpressure is not sufficiently long, the SLS 10 17 Pa s viscosity model appears to be approaching an asymptote. Following the removal of the overpressure at t = 7 years, and hence the stress, we observe strain recovery behavior. As both the Maxwell and SLS models contain an elastic strain component, they both fully recover the elastic strain (ε E ) that was incurred at the onset of pressurization, irrespective of the level of deformation reached during the creep period. However, an important distinction between the Maxwell and SLS models is that the Maxwell rheology does not account for time-dependent viscous strain recovery. This results in the ground surface being permanently uplifted relative to the initial conditions and is dependent on both the viscosity and timescale of the creep behavior. Viscous strain recovery allows the SLS models to return to an undeformed state, if given a long enough time period. The recovery behavior of the KV rheology, like the creep behavior, does not have an elastic component, and so we only observe time-dependent viscous strain recovery.

Expansion
For the expansion models, there is the repeated observation of relaxation behavior following the strain increase at t = 3 years, as seen with the constant expansion models (Figure 4). This includes the relaxation of stress for the Maxwell and SLS models and the convergence of the deformation toward a level that is~0.66 and~0.85 times that of the elastic model, respectively. For each of the rheological models, the reversal of the reservoir expansion at t = 7 years imparts a negative strain on the system that is independent of the model viscosity. The inverse strain incurred in the Maxwell and SLS models is equal in magnitude to the elastic strain (ε E ) incurred at t = 3 years, whereas the response of the KV model is affected by the instantaneous strain approximation mentioned previously. Most notably, all of the models result in apparent subsidence beneath the original ground level. Following this negative strain subsidence response, the Maxwell and SLS models undergo strain relaxation again; however, in this scenario it is restorative and results in uplift that converges with the elastic solution. Dependent on the time interval that is modeled, and the duration of relaxation, a longer-term subsidence signal may remain for models that have a higher viscosity. The KV models cannot undergo stress relaxation, so instead, the deformation time series decay from exaggerated stress states, resulting from an instant change in strain, toward the elastic solution.

Trapezoidal Forcing
The trapezoidal reservoir evolution (Figure 3d) describes a constant-rate injection of magma or exsolution of volatiles, followed by a constant-rate removal of magma or volatiles through intrusion or eruption, or the cooling and contraction of magma. Between these two stages there is a 4-year period in which the reservoir boundary undergoes no change, which may reflect a state of equilibrium. The resultant time series are displayed in Figure 7, with the upper and lower rows being the pressurization and expansion deformation modes, respectively. With this time function, we observe a sequential combination of the behaviors demonstrated by the ramped ( Figure 5) and continuous (Figure 4) time functions, associated with the up-ramp flank from t = 1 year to t = 3 years and the constant deformation mode present between t = 3 years and t = 7 years. Following this plateau, we see compound recovery and relaxation behaviors with the down-ramp flank, from t = 7 years to t = 9 years.

Pressurization
For the pressurization deformation mode, we observe the same creep superposition as for the ramped reservoir evolution (Figure 5), with the up-ramp flank of the trapezoid over the time period t = 1 year to t = 3 years, which is followed by creep behavior for the duration of the constant overpressure. The downramp flank, over the time period t = 7 years to t = 9 years, results in the superposition of strain recovery and produces subsidence across the models due to the decreasing pressure with each time step without 10.1029/2019JB017832 Journal of Geophysical Research: Solid Earth decaying below the original ground level. Unlike the rectangular reservoir evolution (Figure 6), in which we see each of the Maxwell and SLS models incur and recover the same amount of deformation as the elastic model at the onset and removal of pressurization, we observe viscosity-dependent deformation for each of the rheologies across the flanks of the trapezoid function, caused by the ramped, rather than instantaneous, pressure change. It is evident, with the application of the constant overpressure, that the SLS rheology still asymptotes to a level of deformation that is~1.7 times that of the elastic model, apparently irrespective of the prior reservoir evolution. This demonstrates an upper bound to the degree of deformation experienced by this rheology, which may provide simple conversions between elastic and SLS viscoelastic models in order to constrain source parameters. Additionally, for the SLS models, the evolution of deformation from t = 7 years onward is the direct reversal of what is observed prior, as the strain recovery behavior is analogous to the creep behavior in its formulation. This leads to the inference that a symmetrical reservoir evolution results in deformation that is almost fully reversible, if given a sufficient timescale for the viscous strain recovery. Following the return to a zero overpressure, the Maxwell models remain permanently uplifted. A consequence of the KV creep behavior lacking an instantaneous elastic strain response is that the models exhibit an apparent lag between the reservoir pressure evolution and the resultant deformation, which increases with viscosity.

Expansion
Following the removal of strain at t = 7 years, all of the expansion models demonstrate a reversal of the deformation attained prior to this, a result of the deformation being solely driven by the relaxation behavior. During the period of constant strain, t = 3 years to t = 7 years, we observe the convergence of the Maxwell and SLS time series toward a level that is~0.66 and~0.85 times that of the elastic solution, respectively, as seen in the constant ( Figure 4) and rectangular ( Figure 6) reservoir evolutions. Nonintuitively, the Maxwell and SLS models reveal a greater degree of relaxation-related subsidence for the intermediate viscosity (10 18 Pa s) model over the period of constant strain (t = 3 years to t = 7 years) relative to the low-viscosity (10 17 Pa s) model. This results from a compromise between the greater degree of deformation attained by a higher-viscosity model prior to the constant strain and the associated increase in the relaxation time constant. We see that while the high viscosity (10 19 Pa s) model provides the greatest deformation prior  Figure 4) reservoir functions, followed by a ramped reversal. The deformation time series demonstrate fundamental differences between the stress-based (ΔP) and strain-based (Δα) deformation modes for each viscoelastic rheology, by comparing the top row with the bottom row. Comparing the responses to the same deformation mode (i.e., ΔP or Δα, across a row), for the Maxwell, Kelvin-Voigt, and Standard Linear Solid models, highlight the importance of carefully considering a chosen rheology. The maximum elastic vertical deformation is 3.76 and 3.11 cm for the pressurization and expansion deformation modes, respectively. to the constant strain, the timescale for stress dissipation is too long to produce considerable subsidence. Whereas the timescale of the low-viscosity (10 17 Pa s) model is sufficient for full relaxation, the prior deformation is not great enough. The intermediate viscosity (10 18 Pa s) model provides a combination of prior deformation and rate of relaxation that results in the greatest subsidence over the period of constant strain. As expected, the higher-viscosity models demonstrate the least amount of subsidence due to the increased timescales of relaxation with viscosity. As with the ramped ( Figure 5) and rectangular ( Figure 6) reservoir evolutions, we observe exaggerated deformation from the highest-viscosity KV model, resulting from the enhanced stress state imposed by near-instantaneous strain changes, followed by the decay to the elastic solution.

Discussion
With the increasing uptake of viscoelastic rheologies in the modeling of volcanic deformation, there are numerous examples detailing the use of Maxwell (e.g., Newman et al., 2006;Trasatti et al., 2003; and Standard Linear Solid (SLS) configurations (e.g., Hickey et al., 2016;Le Mével et al., 2016), with the Maxwell rheology proving more popular due to its simple formulation. These studies are primarily focused on episodes of uplift or long deformation time series, and often discern reduced overpressure requirements relative to the corresponding elastic models (e.g., Del Negro et al., 2009;Masterlark et al., 2010;Newman et al., 2006) due to the creep behavior. In contrast, viscoelastic models that utilize a volume change at depth are expected to result in reduced levels of deformation, due to the relaxation of crustal stresses (e.g., . We provide a comprehensive comparison of different viscoelastic models under varied conditions, as several important considerations remain unaddressed, such as the unbounded creep behavior of the Maxwell model and the implications of the stress relaxation and strain recovery behaviors. For example, Currenti (2018) recently compared the influence of the Maxwell and SLS models, but the study was limited to a constant overpressure scenario.

Pressurization Versus Expansion
Our results demonstrate the importance of creep versus relaxation behaviors when comparing popular reservoir deformation modes; consequently, the influence of a viscoelastic rheology is dependent on the way in which the deformation is modeled. Strain-based models, invoking a volume change at depth, are expected to result in reduced deformation relative to the elastic model when incorporating a viscoelastic regime, through the relaxation of crustal stresses (e.g., . The time-dependent nature of this relaxation may result in episodes of apparent inflation-induced subsidence to occur (e.g., Figure 4) and may be utilized to explain deflation-related ground deformation patterns (e.g., . Contrastingly, viscoelasticity in stress-based models, through changes in overpressure, results in amplified deformation due to the time-dependent strain increase imparted by the creep behavior. In turn, this reduces overpressure requirements relative to the corresponding elastic model (e.g., Del Negro et al., 2009;Masterlark et al., 2010;Newman et al., 2006). Further to this, the removal of pressure prompts the recovery of strain, which can result in either permanent deformation or a return to the original ground level, dependent on the rheology used. Due to the ambiguity surrounding the mechanisms of reservoir deformation, it is important to consider how the results of geodetic modeling vary between stress-based (ΔP) and strain-based (ΔV) perspectives (Figures 4-7).
The deformation modes in this study were chosen to be analogous, by calculating the volume change owing to the overpressure within an elastic model. Despite the agreement between these parameters, the elastic models resulted in vertical deformations of 3.76 and 3.11 cm, a discrepancy of~21%, for the pressurization and expansion conditions, respectively. This demonstrates that even in a homogeneous elastic model space, reservoir pressurization results in the concentration of the deformation field. This is due to the imposed stresses 'feeling' the free surface, as shown in Supporting Information Figure S6, and consequently, the deformation field is preferentially distributed above the reservoir. As a result, in order for the radially uniform expansion of a reservoir to produce the equivalent level of vertical deformation, an~21% greater volume change than the 'pressure-equivalent' is required in this model setup. This suggests that the volume change of deformation sources from strain-based (ΔV) elastic models may be overestimated with respect to stress-based (ΔP) elastic models. As a consequence of the free surface, this effect is expected to reduce with increasing depth of the deformation source.
While providing an understanding of how viscoelastic time series are affected by different viscosities, our work presents a simplification by utilizing a homogenous model space viscosity. Heterogeneous crustal and thermal properties in the vicinity of a magmatic system are known to partition the resulting deformation field (Del Negro et al., 2009;Gottsmann & Odbert, 2014;Gregg et al., 2012;Hickey et al., 2015Hickey et al., , 2016, and so a natural continuation of this work is to consider the comparison of viscoelastic rheologies and deformation modes in models that more closely represent reality. This can be carried out with models specific to volcanic systems containing heterogeneous crustal properties, provided by seismic tomography, spatially variable rheological effects through the use of viscoelastic shells within an elastic medium (e.g., Currenti, 2018;Currenti & Williams, 2014;Delgado et al., 2018;Newman et al., 2001;Segall, 2016), or a temperaturedependent viscosity distribution that accounts for crustal geotherms and the perturbation owing to a modeled magmatic source (Del Negro et al., 2009;Gottsmann et al., 2017;Gottsmann & Odbert, 2014;Gregg et al., 2013;Hickey et al., 2016). Further to this, an important consideration is the influence of regional stresses and strain fields (Costa et al., 2011;Currenti & Williams, 2014), such as active extension within the Taupo Volcanic Zone (e.g., Cabaniss et al., 2018), on the observed viscoelastic behaviors and the resultant deformation time series.

Asymptotic Behavior and Pressure Scaling
In the time-dependent reservoir evolutions containing a constant overpressure component (Figures 4, 6, and 7), we observe that the deformation time series for the SLS models converge toward a value that is 1.7 times that of the elastic deformation, for a Maxwell arm weighting (μ 1 ) of 0.5, as also recognized by Del Negro et al. (2009). The magnitude of the observed asymptote is governed both by the applied stress and the weight of the elastic arm (μ 0 ). This suggests that for an elastic model with a given overpressure, which is often deemed to be exaggerated, scaling the elastic pressure requirement by the reciprocal of the normalized asymptote constrains a first-order lower bound ("viscoelastic-equivalent") overpressure estimate. We explore this concept in Figure 8, first by examining how the normalized asymptote is affected by the arm weightings (Figure 8a), followed by the influence of the Young's modulus, viscosity, and source overpressure (Figures 8b-8d, respectively). As the weight of the viscoelastic arm (μ 1 ) tends to 1, we observe amplification of both the magnitude and the time constant of the asymptote (Figure 8a). The timescale of convergence is reduced for a greater Young's modulus ( Figure 8b) and increases with a higher value of viscosity (Figure 8c). The magnitude of the asymptote, however, is unaffected by these parameters. In contrast, changes to the source overpressure do not result in any changes to the normalized asymptote, as seen in Figure 8d. Further to this, the magnitudes of the asymptotes are evaluated for incremental viscoelastic arm weightings (Figure 8e), from which an "overpressure factor" is determined (Figure 8f) by taking the reciprocal. In Figure 8g, we demonstrate that an elastic pressure requirement can be scaled by the overpressure factor, dependent on the relative arm weightings, to produce the same long-term deformation signal. The time constant for the asymptotic creep behavior is a function of the elastic parameters and viscosity of the model space, as well as the product of the relative weightings of the spring (μ 0 ) and Maxwell (μ 1 ) arms. In this scenario, as only the arm weights are varied, the time constant is controlled by the product of the arm weights (μ 0 μ 1 ). Figure 8g illustrates that the time series pertaining to models of variable arm weights, given the same product μ 0 μ 1 (i.e., the μ 1 = 0.4 and μ 1 = 0.6 or the μ 1 = 0.2 and μ 1 = 0.8 models), converge to the elastic model over the same timescale, despite exhibiting different ratios of instantaneous (elastic) and time-dependent (viscous) deformation. From the results presented in Figures 8a-8d, with the magnitude of the asymptote affected only by the arm weights, this overpressure scaling and determination of a viscoelastic-equivalent pressure requirement is universally applicable. While the models presented within this investigation contain major simplifications when compared to volcanic settings (i.e., homogeneous crustal parameters), they provide, as intended, insight into the expected rheological response to a deformation episode, prior to incorporating increasing levels of complexity such as spatially variable crustal properties. Variations in viscosity and mechanical properties within the vicinity of a magmatic system, owing to the local thermal regime and the extent of crustal structure and heterogeneity, can greatly affect the partitioning of strain (Del Negro et al., 2009;Gottsmann & Odbert, 2014;Gregg et al., 2012;Hickey et al., 2015Hickey et al., , 2016, resulting in spatial patterns and timescales of deformation that are individual to each volcanic center for a given change in reservoir conditions.

Modeling of Subsidence Episodes
In the above results and analysis, we focus on the application of an overpressure or expansion, in which the resultant stresses and strains are directed outward from the centroid of the modeled reservoir. Conversely, (e) Normalized asymptote (NA) for incremental weightings of the Maxwell arm (μ 1 ). (f) The value of the normalized asymptote (NA; from (e)) can be used to calculate an overpressure scaling factor (OF) for different arm weightings, where OF = 1/NA. (g) The scaling of overpressure requirements for any given arm weightings allows first-order viscoelastic constraints to be determined from a given elastic model. Overall, the weight of the elastic (μ 0 ) arm determines the magnitude of the instantaneous elastic response, and so the proportion of the viscous creep component increases with the increasing weight of the Maxwell (μ 0 ) arm. The time constant of the creep behavior, in this parameter space, is controlled by the product of the arm weights (μ 0 μ 1 ). directing the stresses and strains toward the centroid of the modeled reservoir, underpressure or contraction conditions are produced. Due to the nature of the time-dependent behaviors (section 2), reversing the direction (i.e., changing the sign) of the stress or strain condition will result in the same behaviors, but with a negative magnitude (as seen in Supporting Information Figures S2-S5). Normalizing these subsidence time series with respect to the corresponding elastic models produces the same plots as shown in Figures 4-7. Consequently, our observations relating to uplift can be extended to a subsidence scenario, and so we recognize that viscoelasticity still reduces the pressure requirement of an elastic model. However, the loss of volume at depth can result in apparent ground subsidence induced by source inflation, due to the relaxation of crustal stresses.

Viscoelastic Behaviors
The prior analysis of the constitutive viscoelastic behaviors, coupled with our model results, demonstrates that there are numerous behavioral characteristics belonging to the Kelvin-Voigt (KV) rheology that make it unsuitable for the modeling of crustal materials and volcano deformation. These include the lack of an elastic creep response and being unable to accommodate stress relaxation, as well as the infinite stress response to a change in strain. Compared to the other configurations, KV viscoelasticity is rarely used. In contrast, the Maxwell rheology is widely utilized in pressurization models due to its simple formulation, but only provides a first-order representation of the creep and strain recovery behavior. The linearly unbounded creep response, combined with the lack of time-dependent strain recovery, ultimately allows large irreversible strains to be generated. This is true even for small stresses, if they are applied for a sufficiently long time period. The SLS rheological model, however, is the simplest viscoelastic configuration to provide a reasonable representation of the time-dependent behaviors. This includes creep behavior that features both an instantaneous (elastic) and an asymptotic time-dependent (viscous) component, an asymptotic stress relaxation response, and the ability to fully recover induced strains over a sufficiently long time period. With the SLS model, we have demonstrated that the relative weightings of the elastic arm (μ 0 ) and the Maxwell arm (μ 1 ) strongly influence both the asymptote magnitude and the timescale of convergence; however, there is little experimental evidence to suggest how these values should be allocated other than on the basis to reduce misfit. As a result, an equal weighting of 0.5 is most common (e.g., Del Negro et al., 2009;Gottsmann & Odbert, 2014;Hickey et al., 2013Hickey et al., , 2015Hickey et al., , 2016Le Mével et al., 2016;Morales Rivera et al., 2018). However, insight into the weighting of the fractional shear moduli may be gained from the ratio of instantaneous (elastic) to time-dependent (viscous) deformation (Figure 8g), and so it may be possible to derive alternative relative contributions from deformation time series, but the use should be carefully considered.

Implications for the Interpretation of Monitoring Data
While this study focuses on forward modeling and the resultant time series, many investigations determine the parameters of deformation sources, including their pressure and/or expansion histories, through the inversion of geodetic data. Based on the differences between creep behaviors observed in the above results, we consider how the choice of the implemented rheology may affect the inference of the pressure evolution of a magmatic system. In Figure 9, we consider the interpretation of a linearly increasing deformation time series based on the implemented rheology, for a given set of model parameters. If modeled elastically, the pressure evolution directly reflects the deformation data, and so a steadily pressurizing system is inferred. However, a corresponding model using the Maxwell rheology may suggest that the system is instead maintaining its overpressure, due to its linear creep behavior, whereas the asymptotic nature of the SLS creep behavior may suggest that the rate of pressurization is decreasing with time. It is clear, from this hypothetical situation, that the response to an episode of unrest is critically dependent on the implemented rheology, foremost with the Maxwell model suggesting that the system is undergoing no change and the elastic model implying an amplified rate of pressurization. Determining the temporal evolution of a reservoir during episodes of unrest remains an ever-present challenge for the monitoring of volcanic systems, due to complexities within observed deformation time series, heterogeneous crustal properties and structures, and rheological assumptions. Despite often satisfying a goodness-of-fit criterion between observed ground deformation patterns and model predictions, the applicability of purely elastic behavior in the vicinity of a subvolcanic system is likely limited. The thermal evolution of long-lived magmatic systems (e.g., Annen, 2011;Gelman et al., 2013;Karakas et al., 2017) is more consistent with a nonelastic middle and upper crustal rheology, containing a time-dependent (viscous) component of behavior, and may better represent the observed deformation field (e.g., Jellinek & DePaolo, 2003;Newman et al., 2001). Nonlinear relationships between reservoir evolutions and the observed deformation time series in viscoelastic models necessitate the construction of models that closely represent individual volcanic systems, through the inclusion of heterogeneous thermomechanical crustal properties, to better understand the importance of timedependent rheological behaviors when modeling volcanic ground deformation. As the exact rheological behavior of crustal materials in volcanic settings remain ambiguous, it is important to consider the evolution of a magmatic system from multiple viewpoints.

Conclusions
In this study, we have analyzed the fundamental, time-dependent behaviors belonging to the Maxwell, Kelvin-Voigt (KV), and Standard Linear Solid (SLS) viscoelastic configurations, and demonstrate that the influence of the rheology is dependent on both the modeled reservoir evolution and the implemented deformation mode. Stress-based, or pressurization, models generally result in amplified deformation with the inclusion of a viscoelastic rheology, which in turn allows the pressure requirements of an elastic model to be reduced. Strain-based, or volume change, models that are commonly produced through analytical inversions will experience reduced deformation relative to the elastic result. These models may also exhibit episodes of apparent inflation-induced subsidence or subsidence-induced inflation, related to the dissipation of crustal stresses. While the inferred source evolution from classical elastic deformation models directly reflects the profile of a deformation time series, we demonstrate that time-dependent viscoelastic behaviors can produce deformation time series that deviate significantly from the profile of the modeled reservoir evolution. Consequently, we establish that the characterization of an unrest episode is critically dependent on the rheology utilized in the deformation model and further demonstrate that determining the evolution of a reservoir from observed deformation patterns remains an ever-present challenge in geodetic modeling. Moreover, by analyzing the behaviors exhibited by the commonly used Maxwell rheology, we suggest that despite its simple formulation, it is largely unsuitable for the modeling of volcanic deformation due to the capacity for generating large irreversible strains. While the SLS rheology is believed to most accurately represent the time-dependent behaviors of crustal materials, and perhaps offers the best combination of simplicity and realism, the relative contribution of the elastic and Maxwell arms is an issue to consider. Figure 9. Schematic demonstrating the variation in the inferred pressure evolution of a magma reservoir, based on the implemented rheology for a given model setup. The shaded area represents the solution space for the Standard Linear Solid model based on the elastic and Maxwell end-members, dependent on the weighting of the Maxwell arm (μ 1 ).