Time-dependent ground deformation is a key observable in active magmatic systems, but is challenging to characterize. Here we present a numerical framework for modeling transient deformation and stress around a subsurface, spheroidal pressurized magma reservoir within a viscoelastic half-space with variable material coefficients, utilizing a high-order finite-element method and explicit time-stepping. We derive numerically stable time steps and verify convergence. We then explore the frequency dependence of surface displacement associated with cyclic pressure applied to a spherical reservoir beneath a stress-free surface. We consider a Maxwell rheology and a steady geothermal gradient, which gives rise to spatially variable viscoelastic material properties. The temporal response of the system is characterized numerically with a transfer function that connects peak surface deformation amplitude and phase lag with respect to sinusoidal reservoir pressurization. Amplitude and phase of this transfer function vary with the frequency of pressure forcing. The volume of host rock exhibiting viscous response around the reservoir is also frequency dependent, depending on a threshold local Deborah number that measures the characteristic timescale for pressurization against a spatially varying Maxwell relaxation time. Although commonly idealized as a thin shell, we find that this dominantly viscous region may define a spatially complex brittle-ductile transition, depending on crustal thermal state, at longer forcing periods. Because arbitrary chamber pressure histories can be represented through a superposition of sinusoidal forcing, our results and general computational framework apply to a wide range of transient deformation scenarios relevant for characterizing transcrustal magmatic systems.

Time‐dependent ground deformation is a key observable in active magmatic systems, but is challenging to characterize. Here we present a numerical framework for modeling transient deformation and stress around a subsurface, spheroidal pressurized magma reservoir within a viscoelastic half‐space with variable material coefficients, utilizing a high‐order finite‐element method and explicit time‐stepping. We derive numerically stable time steps and verify convergence, then explore the frequency dependence of surface displacement associated with cyclic pressure applied to a spherical reservoir beneath a stress‐free surface. We consider a Maxwell rheology and a steady geothermal gradient, which gives rise to spatially variable viscoelastic material properties. The temporal response of the system is quantified with a transfer function that connects peak surface deformation to reservoir pressurization in the frequency domain. The amplitude and phase of this transfer function characterize the viscoelastic response of the system, and imply a framework for characterizing general deformation time series through superposition. Transfer function components vary with the frequency of pressure forcing and are modulated strongly by the background temperature field. The dominantly viscous region around the reservoir is also frequency dependent, through a local Deborah number that measures pressurization period against a spatially varying Maxwell relaxation time. This near‐reservoir region defines a spatially complex viscous/elastic transition whose volume depends on the frequency of forcing. Our computational and transfer function analysis framework represents a general approach for studying transient viscoelastic crustal responses to magmatic forcing through spectral decomposition of deformation time series, such as long‐duration geodetic observations.