History-Dependent Volcanic Ground Deformation From Broad-Spectrum Viscoelastic Rheology Around Magma Reservoirs

of

as the scales of deformation vary.A spectral domain perspective on unsteady magmatism has been applied to eruption sequences (Jellinek et al., 2004), but has yet to be explored in the context of active geodetic observations at volcanoes.Under an appropriate description of magma chamber deformation, for example, physical processes become linear or nonlinear spectral "filters" which modify an input signal (e.g., a magmatic intrusion) to generate frequency and wavenumber spectra of output signals (e.g., ground displacement).
We derive an analytic transfer function between magma chamber pressure and maximum surface deformation, extending the findings of Rucker et al. (2022) that frequency-dependent phase and amplitude are signatures of the viscoelastic filter.These metrics are shown to be sensitive to magma chamber geometry, host rock rheology, and crustal temperature field.This approach establishes a surprisingly simple link between analytic (e.g., Dragoni & Magnanensi, 1989;Segall, 2016) andnumerical (e.g., Del Negro et al., 2009;Gregg et al., 2013;Head et al., 2021) viscoelastic magma chamber models.
We then demonstrate through two volcanologically-motivated examples that viscoelastic host rock response to broadband magma chamber pressurization-by which we mean a magma pressure time-history whose power spectrum contains many frequency components-leads to history-dependence of surface deformation.This has implications for the interpretation of geodetic timeseries, inference of crustal stress state, and eruption forecasting.We show that the elastic/viscous transition point around magma chambers, a key indicator of magma storage maturity (Degruyter & Huber, 2014), arises via frequency-dependent partitioning of viscous and elastic crustal strains when crustal rheology varies with temperature.

Frequency Domain Viscoelastic Model for Magma Reservoir Deformation
We consider a spherical magma chamber with radius r o and depth d in a viscoelastic half space (Figure 1a).The magma chamber has a uniform interior temperature T in .Steady state crustal temperature T(r) is assumed radially symmetric (Figure 1a), determined by two boundary conditions T(r o ) = T in and T(d) = 0.This temperature field approximates a crustal geotherm dominated by the thermal perturbation of a magma reservoir, and becomes inaccurate at distances where the radial gradient dT/dr approaches the background vertical gradient.For a magma chamber with interior temperature of 1,000-1,200°C at depth of 5-10 km and vertical geothermal gradient of 20-30°C/km, a radially symmetric temperature profile approximates an extended region above the chamber adequately to a radius equal to the depth for some cases (Figure S1 in Supporting Information S1).
Rheology of the crust is assumed Maxwell viscoelastic, with temperature-dependent viscosity that increases exponentially away from the chamber (Figure 1a and in Supporting Information S1) (Del Negro et al., 2009;Karlstrom et al., 2010).This rheological model has been used extensively in volcano and tectonic geodesy for understanding viscoelastic response (e.g., Segall, 2010;Tromp & Mitrovica, 1999).However it is not unique and likely over-simplifies the parameterization of microscale creep mechanisms that set the relative proportions of recoverable versus non-recoverable strains for given forcing (H.C. Lau & Holtzman, 2019).We use the Maxwell model for demonstrating a theoretical framework that is transparent and extendable.
Mechanical response to forcing associated with chamber pressurization, such as surface ground deformation, is governed by the rheology of the magma reservoir plus host rock system.Time-dependent pressure forcing is a boundary condition for radial stress at r o , which, together with an implicit boundary condition of vanishing displacement at r = ∞, predict a unique displacement pattern.Surface displacement is found via an approximate correction on the full space solution for stress-free surface conditions (Del Negro et al., 2009;Mogi, 1958).
It is often assumed that the geometry of magma reservoirs is fixed in time and that magmatic forcing arises from magma pressurization relative to lithostatic stress.With these assumptions, viscoelastic deformation around a magma reservoir is a Linear Time-Invariant system (Schetzen, 2003) and can be solved via transfer function in the frequency domain between the forcing function as input signal and any resultant scalar mechanical output such as vertical displacement at a point on the Earth's surface.Although in general this transfer function can be found numerically (Rucker et al., 2022), we develop an analytic transfer function for the model problem above (derived in Supporting Information S1) to demonstrate key phenomenological outcomes of broad-spectrum viscoelastic deformation.
The transfer function  {()|()} = f∕ f linearly relates the frequency spectrum of an unknown quantity  f (i.e., an output signal, where hat signifies Fourier transform) to an input signal  f .This frequency domain method provides an equivalent approach to time-domain methods for forward and inverse problems (Figure 1b).While any input or output signals can be used to compute   , we will assume for illustrative purposes that the chamber overpressure P(t) relative to lithostatic is a known input signal, and the maximum surface displacement u z (t) is the output signal (Figure 1a).
For a spherical magma chamber in an elastic halfspace the transfer function, assuming free surface corrections (McTigue, 1987;Segall, 2016) appropriate for chambers with r < D/2.5, may be stated as where superscript el refers to linear elastic rheology, μ is shear modulus and K is the bulk modulus.The elastic transfer function is hence real-valued and independent of frequency.We are interested in isolating viscoelastic effects, and henceforth normalize transfer functions by Equation 1.
Applying the constitutive relations, boundary conditions, and Fourier transform (supplement), we obtain the normalized viscoelastic transfer function   between chamber pressure and maximum surface displacement This equation introduces the Deborah number De(r) = ωη(r)/μ, a dimensionless product of characteristic cyclic strain rate (ω > 0) or inverse period (τ = 2π/ω) associated with reservoir pressurization, and Maxwell time η(r)/μ that measures the relative importance of viscous versus elastic response in the system.De varies with space when material properties η and/or μ vary with space, here due to the Arrhenius relationship between viscosity and temperature (Figure 1a).In general, the elastic moduli K and μ also vary with temperature, but for temperature difference between 0° and 800°C K varies by 1% and μ varies by 10% (Bakker et al., 2016).This is small in comparison to the orders of magnitude of variation in viscosity, so we assume uniform elastic bulk moduli K and μ to keep the model simple.
From Equation 2we see that the normalized transfer function   is a viscoelastic correction on the elastic displacement.The magnitude of this correction is determined by the spatial structure of Deborah number.If the crust is elastic, the viscosity is infinitely large and De(r) = 0 for all frequencies; in this case the normalized transfer function becomes unity.
Because crustal magmatic systems are characterized by spatially localized temperature anomalies, the transfer function defines a localized region of viscous response and implies that a near-chamber elastic/viscous transition controls viscoelastic deformation.This is the motivation behind magma chamber models that assume a discrete viscoelastic shell of uniform viscosity (Degruyter & Huber, 2014;Dragoni & Magnanensi, 1989;Karlstrom et al., 2010), a simplification of continuous variation in material properties expected in spatially variable thermal field.
We find that the general response of a material with variable viscoelastic moduli can be precisely represented by that of a discrete and uniform viscoelastic shell with frequency dependent effective viscosity and temperature under monochromatic forcing at period τ.For such case, the outer radius of the viscoelastic shell R eff defining a transition to elastic response, the uniform effective temperature T eff and the resulting (constant) Deborah number De eff associated with viscosity η eff (T e ff) are found by requiring that the transfer function matches Equation 2. Generic unsteady pressure forcing can be decomposed into superpositions of such harmonic forcing functions (Rucker et al., 2022).This result thus provides a bridge between classical analytic models and numerical approaches for magma chamber deformation.

Properties of the Transfer Function
Polar decomposition of the complex-valued   calculated by Equation 2via   =   allows for the identification of an amplification factor   = || and a phase delay of −φ.These quantities completely characterize the transfer function, and provide insight into the frequency dependence of viscoelastic deformation.
Both the amplification factor and the phase delay frequency spectra vary with the physical properties of the system.Figures 2a and 2b shows dependence on the two parameters controlling steady-state temperature of the crust, the chamber temperature and assumed background thermal gradient (defined using a radius contour associated with chamber depth below surface).As shown in Figure 2b,   increases with forcing period due to increasing spatial extent of viscous response.In contrast, the frequency dependence of the phase lag is non-monotonic (Figure 2a), and peaks at a value well below the Newtonian fluid phase lag of π/2 (Figure 2d).Time-domain delay between deformation and harmonic pressure forcing increases monotonically with period (Figure S5 in Supporting Information S1).
For short forcing periods, the crust becomes more elastic hence the phase lag decreases with forcing period.As forcing period τ = 2π/ω increases, the phase lag reaches a maximum and drops to zero as τ → ∞.This behavior can be explained by a competition between spatially confined viscous relaxation and cyclic pressure forcing.The spatial extent of viscous response R eff is confined near the reservoir by the viscosity profile (set by temperature).As forcing period increases, effective temperature of the shell decreases while R eff increases (Figure 2c).The timescale for viscoelastic stress relaxation over the shell   ∕( ∕0) 3 (Dragoni & Magnanensi, 1989) thus increases with forcing period, although shell effective Deborah number is dominated by τ and decreases monotonically.The increase of   with τ reflects viscous strain across the growing shell, with maximum viscoelastic response implicated by a maximum phase lag φ.For fixed R eff , maximum phase lag occurs where   =  ( ∕0) −3∕2 (Figure S6 in Supporting Information S1).Large forcing period becomes an effective elastic system (zero phase lag) with reservoir radius R eff rather than r 0 (Karlstrom et al., 2010), within which viscous stresses are uniform.
The spatial viscosity structure that gives rise to this behavior, perhaps a defining feature of magma reservoirs, means that Newtonian viscous crustal response is never realized in our model (red lines in Figure 2d).As shown in the Supplement, if host rock viscosity is constant,   does approach the Newtonian limit at large τ.Of course, at long forcing periods and spatial scales much larger that the reservoir we expect that the vertical variation in temperature controls viscous response (e.g., Tromp & Mitrovica, 1999) and our model assumptions are no longer valid.Additional structure in the phase lag spectra may also arise from details of the problem not modeled here (Rucker et al., 2022), and our results are specific to Maxwell rheology.But our analytic model captures first-order frequency dependence of the reservoir-crust system on timescales of interest to volcano geodesy.
We also see that crustal thermal structure is imprinted on the transfer function.Chamber depth is approximated by the distance d between surface temperature and chamber boundary T in .Larger T in and d generally result in more pronounced phase lag, with chamber temperature being particularly sensitive (Figure 2a, Figure S4 in Supporting Information S1), and noting trade-offs in maximum phase lag between these parameters.Longer forcing periods (lower frequency) require larger effective shells with smaller effective temperatures.Likewise, a hotter crust requires a larger effective shell with larger effective temperature.Although not explored here explicitly, the trajectories in Figure 2d show how variable material coefficients and constitutive model interact: in this case, the background thermal field dictates the range of possible Maxwell viscoelastic response for harmonic forcing functions.
Of course, the discrete shell equivalence explored in Figure 2 is only precise when the pressure forcing is sinusoidal in time.For broadband pressure forcing with more than one frequency as will be explored in the next section, there does not exist a set of values (R eff , T eff ) that will allow the effective shell to yield the same response for all frequency components.

Signatures of Viscoelastic Volcano Ground Deformation
It is useful to demonstrate the utility of the spectral transfer function approach through application to synthetic time-series that mimic geodetic observations.In this section we show examples of maximum vertical surface displacement and crustal strains obtained using the spectral method for various pressurization episodes of a subsurface reservoir.

Isolated Pressurization Events
We consider two idealized forcing functions, a square pulse and sawtooth pressurization history of the chamber above some baseline.Sawtooth pressurization approximates gradual inflation and rapid deflation inferred at many erupting volcanoes (e.g., Chadwick Jr et al., 2022;Delgado et al., 2018;M. R. Patrick et al., 2020), while a square pulse provides a contrasting impulsive forcing with greater long-period energy.Elastic surface displacement would be proportional to this chamber pressure (Segall, 2010), however viscous creep causes the displacement to be larger and last longer (Figure 3a).These differences are evident in the displacement frequency spectra (Figure 3a inset), as the viscoelastic transfer function increasingly amplifies periods longer than the characteristic time of the pulse.Such frequency-dependent amplification results in larger displacement and continuing surface deformation after the end of the episode, a larger effect for the square pulse due to increased energy at long periods relative to the sawtooth pulse.Despite the broadband forcing function, net viscoelastic response in surface displacement thus varies with spectral content of the pressurization.As shown in Figure S3 in Supporting Information S1, a rapid pressure forcing event (e.g., square or sawtooth episode with duration of 1 day) exhibits surface displacement close to the elastic solution.Conversely for longer duration forcing, displacements are amplified and the relaxation lags input pressure.
Similar to the harmonic forcing case, the physical origin of this mechanical response can be understood by examining viscous strain around the chamber and the location of effective elastic transition.For linear viscoelastic rheology, the deviatoric strain is the sum of an elastic and a viscous component.We calculate the partition of the total viscous deviatoric strain in space     ∕   , which decreases away from the chamber with the background temperature of the crust (Figure 3b colored lines).Viscous relaxation is thus focused around the chamber, with effectively elastic deformation far from the chamber.We characterize the point where     ∕   < 10% as an elastic/viscous transition point for the variable coefficient mechanical response to chamber pressurization.With different magmatic parameters and pulse durations,     ∕   appears to vary spatially (Figure S7 in Supporting Information S1).However, the partition between the viscous and elastic components collapses onto a single curve when plotted against local Deborah number (Figure 3b).This data collapse is a manifestation of the "thermorheologically simple" viscoelastic behavior (Muki & Sternberg, 1961) of our constitutive model.Rucker et al. (2022) suggest that the extent of significant viscous strain is characterized by a contour where De(r) ≈ 10, for harmonic forcing.Figure 3b shows that this result does not extent to broadband forcing, with the sawtooth transition point at De(r) ≈ 35 and square pulse at De(r) ≈ 65.The frequency spectrum of pressurization is thus imprinted on the spatial pattern of viscoelastic response.

The dependence of transfer function 𝐴𝐴
 on control parameters suggest possible constraints on thermorheologic and geometrical aspects of the combined magma reservoir-crust system (defining the top of a transcrustal magma transport network, Sparks et al. ( 2017)) using frequency domain observations.With background temperature as an example, we see in Figure 2 that the 1,200°C chamber always amplifies the elastic displacement more strongly than the 800°C chamber.However, the larger relative phase delay between models switches as forcing period increases.The maximum difference in phase delay between 1,200°C chamber and 800°C chamber occurs around 3.8 days (Figure S4 in Supporting Information S1).For broadband forcing functions, such signatures imply specific relationships between deformation frequencies.

History Dependence in Event Sequences
The frequency dependence of viscoelastic response also implies that sequences of impulsive pressure forcing functions with long period energy may exhibit incomplete viscous stress relaxation that imparts history dependence to the resulting deformation.We demonstrate this time-dependent stress build-up by considering ground deformation response to sawtooth pulses of duration t p in sequence (Figure 4).In the case shown in (Figure 4a), the repose time t s (e.g., inter-event or hiatus time) between each pressure pulse increases through the sequence.The resulting surface uplift consists of repeated episodes, each with a different peak vertical uplift, despite identical forcing functions.Larger peak displacements occur when the repose time t s is less than the pressurization pulse duration t p .But even for the last event in the sequence, occurring after a repose time of twice the pulse duration (100 days), there is amplification relative to an isolated pulse (the amplitude of the first uplift event).
We explore this relationship in more detail in Figure 4b using two sawtooth pressure pulses, each with similar t p but variable t s .Peak ground deformation relative to an isolated pulse is greater than 1% even for t s /t p of 30 (e.g., representing ∼3 months episodes separated by ∼8 years).This history dependent behavior is connected again to the extent of viscous relaxation in the crust at the onset of each pressure forcing episode, which involves both the pressurization and repose history.For short repose time, the crust does not have sufficient time to relax before subsequent pressurization occurs, hence the pre-stressed crust results in larger displacement at the end of the second pressurization episode.For large t s , crustal stress relaxes over a greater range of periods, hence the peak displacement is closer to the previous episode.For the two square pulse case (Figure 4b), the peak surface uplift in the second episode is always higher than the peak uplift during the first episode, which corresponds to a zero-stress, fully-relaxed initial condition.The range of periods for which relaxation can occur is set by the spectrum of Maxwell relaxation times in the crust, determined by the thermal profile in our model.The product of forcing function by the transfer function   at each period determines its relaxation time.Forcing function with energy in longer period and hotter crust result in larger amplitude history dependence of surface deformation (Figure S8 in Supporting Information S1).
Incomplete relaxation of forcing signals with significant long-period power is thus likely on timescales of repose times at active volcanoes (e.g., years to centuries).This history dependence implies that the absolute stress state of the crust-a key initial condition for chamber failure or eruption triggering-reflects concurrent magmatism but also magmatic and tectonic history.

Summary and Implications
Viscoelastic deformation has been long recognized as a likely contributor to volcano ground deformation.The spectral method introduced here is a powerful framework for studying time-dependent deformation in magmatic systems where linear viscoelastic crustal rheology can be assumed.In such systems, because each frequency component evolves independently, a general surface deformation signal may be decomposed into constituent frequency spectra that we have shown vary with crustal temperature and geometry of storage.While this approach does not resolve non-uniqueness in ground deformation, it suggests new analysis tools for studying restless volcanoes.In particular, we have shown that viscoelasticity leads to a frequency-dependent lag between source pressure and surface deformation.Stress state and surface deformation associated with intrusion sequences is therefore likely history dependent.
Our approach extends analytic techniques for obtaining viscoelastic solutions from corresponding elastic solutions (e.g., Bonafede & Ferrari, 2009;Dragoni & Magnanensi, 1989;Fung, 1965;Segall, 2016) by including variable material properties and an approximate free-surface correction.As shown in Rucker et al. (2022), transfer functions can be obtained numerically, which allows for arbitrary background thermal and mechanical material heterogeneity, chamber geometry, topography/layering, and tectonic stresses.The advantage of our analytic approach is its physical transparency and efficiency of implementation.Except in unusual cases (M.Patrick et al., 2019), pressure itself is an unknown quantity (Anderson et al., 2020;Zhan et al., 2019) so reservoir source-time forcing function as well as the transfer function must be inferred from observations.Multiphase magma mass balance and buoyancy impart additional complexity to viscoelastic deformation in general (Segall, 2019;Sigmundsson et al., 2020).
Of course, quantitative details of our results depend on assumed rheology and thermal structure.Here we use a simple radial temperature profile, which approximates expected steady near-surface gradients (Figure S1 in Supporting Information S1) and localized viscoelastic rheology associated with a thermal boundary layer around the chamber.Numerical explorations (e.g., Gregg et al., 2013;Karakas et al., 2017) show that this solution is a much poorer approximation below the chamber, which will lead to significant deviation of predicted viscoelastic deformation at long periods.Although Rucker et al. (2022) demonstrate that the frequency domain structure of transfer function  {| } found numerically is qualitatively similar to our analytic results, dependence on the spatial distribution of material parameters, domain geometry, and far-field boundary conditions have yet to be established.Additional time-dependent processes, such as poroelasticity (Liao, 2022;Liao et al., 2018;Mittal & Richards, 2019), or other constitutive models (e.g., Head et al., 2021), could be incorporated into the spectral method developed here.
It is also important to evaluate our assumptions of steady state geometry and background temperature.On timescales of sustained magmatism, active volcano lifetimes (100s Kyr for stratovolcanoes, Hildreth, 2007) are smaller than thermal conduction timescale based on chamber depth (d 2 /κ ∼ 800 kyr for d = 5 km and thermal diffusivity κ = 10 −6 m 2 /s), although often longer than a conduction time based on intrusion radius of 1-30 Kyr (for R 0.2-1 km).This range provides an upper bound on variants in the thermal field.Temperature anomalies that precede eruptions (Girona et al., 2021) suggest far shorter timescales of thermal transience mediated by fluid advection.However, if conduction-dominated, temperature may be reasonably approximated as quasistatic (although likely not steady state) on ∼1 kyr timescales or less.Magma reservoir geometry also evolves on a wide range of timescales (e.g., M. R. Patrick et al., 2020;Rivalta et al., 2019), although large changes in magma storage geometry may be linked to thermorheologic evolution of the crust (Colón et al., 2018;Huber et al., 2019;Karlstrom et al., 2017).Such temporal variations violate the Linear Time-Invariant nature of the transfer functions presented here without generalization of the model.For example, evolution of temperature or other scalar fields in the system would require spectral treatment of a more complete set of governing equations.These caveats aside, we believe that the frequency domain framework outlined here has potential both for observational validation of viscoelastic volcano models, and for enhancing understandings of the diverse deformation patterns seen at volcanoes globally.To argue the first point, we note that the most significant variation in the transfer function occurs in a range of periods (1-100s of days) where geodetic observations are common.This leads to a hypothesis that viscoelasticity is a common but under-recognized feature of many volcano geodetic signals.For example, apparent variation in inflation thresholds invoked for eruption cycles at Axial seamount (Chadwick Jr. et al., 2022) could instead reflect imbalanced viscoelastic relaxation of host rocks at long periods relative to short.By leveraging the full displacement field (not just maximum vertical displacement presented here), we expect that viscoelastic transfer functions can be validated with long-term geodetic timeseries at volcanoes.
As for the second point, we note that the global distribution of active volcanoes represent diverse snapshots of slowly evolving transcrustal magmatic systems.Identifying quantitative ways in which a thermally immature versus mature volcanic system should deform as a function of frequency provides concrete predictions that can be used to design experiments and study integrated deformation recorded by ductile strain markers (e.g., McNulty et al., 1996) seen around exhumed magmatic systems.

Figure 1 .
Figure 1.(a) Geometry of the model, representative vertical surface deformation (top black curve), temperature (inset black curve) and viscosity profiles (yellow curve).(b) Representation of the spectral method and transfer function approach connecting general input and output signals f in and f out .The input signal is represented either in time (t) or frequency domain (ω), hats represent Fourier transform   .Multiplication by transfer function  {|} is the frequency domain equivalent of time domain solution to the linear governing equations.

Figure 2 .
Figure 2. (a) Phase lag φ and (b) amplification factor   of the normalized transfer function   calculated according to Equation 2 as functions of harmonic pressure forcing period (τ = 2π/ω), for different chamber temperatures and depths approximated from a thermal model.(c) Effective outer radius R eff and temperature T eff of equivalent viscoelastic shell for harmonic pressure forcing.(d) Phase lag φ of the transfer function as function of effective viscoelastic shell radius and Deborah number in the shell (De eff , colored contours).The red curves are trajectories from variable coefficient models shown in (c), illustrating the dependence of phase lag on shell Deborah number and size.The Newtonian fluid limit (φ = π/2) is not reached by the red lines because viscoelastic relaxation saturates in the variable coefficient system, with effective elastic response observed at small τ and approached as τ → ∞.

Figure 3 .
Figure 3. (a) Maximum vertical ground displacement u z associated with sawtooth and square pulse-shaped pressurization over 100 days with reference viscoelastic parameters (Table S1 in Supporting Information S1) and d/r o = 5, T in = 1,000°C.Dashed lines show the corresponding elastic case.Insert shows amplitude spectra of the normalized displacement obtained from the transfer function as functions of period τ.(b) Spatial extent of viscous relaxation for different pressurization time-histories, measured by the ratio of maximum deviatoric viscous strain     and maximum total deviatoric strain     , plotted against spatially variable Deborah number De(r).Horizontal gray line indicates     ∕   = 10% , an elastic/viscous transition point.These curves are insensitive to chamber temperature and depth (Supporting Information S1), demonstrating that the rheologic transition depends on spectral content of the forcing function: De ≈ 35 for sawtooth forcing (black curve), De ≈ 65 for square pulse (red curve), De ≈ 10 for harmonic pressurization (dotted yellow curve).

Figure 4 .
Figure 4. (a) Maximum vertical ground displacement (solid line, left axis) as a function of time in response to sawtooth pressure forcing (dash line, right axis) and reference model parameters.Each pressure forcing episode has duration t p = 100 days.The durations of repose time between subsequent pulses are t s = 25, 50, 100, and 200 days.(b) Grounddisplacement from two pressurization episodes with variable repose time t s .Each pulse has a duration of t p = 100 days.Insert shows the relative increase of the peak uplift in the second episode with regard to the peak uplift in the first episode and decreases with repose time.The system has reference viscoelastic parameters (TableS1in Supporting Information S1), d/ r o = 5, T in = 1,000°C.