Decompression and Fracturing Caused by Magmatically Induced Thermal Stresses

Studies of host rock deformation around magmatic intrusions usually focus on the development of stresses directly related to the intrusion process. This is done either by considering an inflating region that represents the intruding body, or by considering multiphase deformation. Thermal processes, especially volume changes caused by thermal expansion are typically ignored. We show that thermal stresses around upper crustal magma bodies are likely to be significant and sufficient to create an extensive fracture network around the magma body by brittle yielding. At the same time, cooling induces decompression within the intrusion, which can promote the appearance of a volatile phase. Volatile phases and the development of a fracture network around the inclusion may thus be the processes that control magmatic‐hydrothermal alteration around intrusions. This suggests that thermal stresses likely play an important role in the development of magmatic systems. To quantify the magnitude of thermal stresses around cooling intrusions, we present a fully compressible 2D visco‐elasto‐plastic thermo‐mechanical numerical model. We utilize a finite difference staggered grid discretization and a graphics processing unit based pseudo‐transient solver. First, we present purely thermo‐elastic solutions, then we include the effects of viscous relaxation and plastic yielding. The dominant deformation mechanism in our models is determined in a self‐consistent manner, by taking into account stress, pressure, and temperature conditions. Using experimentally determined flow laws, the resulting thermal stresses can be comparable to or even exceed the confining pressure. This suggests that thermal stresses alone could result in the development of a fracture network around magmatic bodies.


of 25
Studies of host rock deformation around magma chambers usually focus on stresses directly related to magma transport (such as dyke or sill emplacement). The customary approach is to prescribe the magma body as an overor underpressured volume, representing an inflating or deflating region within the crust. There are analytical solutions that describe the displacement or stress field for different intrusion geometries in a purely elastic host rock (Fialko et al., 2001;Kiyoo, 1958;McTigue, 1987;Yang et al., 1988). However, if large volume changes are considered, equivalent of more than several MPa pressure difference, a few km below the surface, brittle failure becomes increasingly likely due to the small confining pressure. In this case, a purely elastic rheology is no longer applicable and the quantification of the stress state and tensile or dilational shear failure is of particular importance. This is because fractures or dikes propagating from the inclusion might reach the surface, resulting in an eruption or in the appearance of fumaroles. To investigate stresses and deformation in a visco-elasto-plastic host rock, several studies applied thermo-mechanical numerical modeling. Some utilize a visco-elastic rheology to quantify stresses and determine the onset of failure (e.g., Gregg et al., 2012;Head et al., 2022;Novoa et al., 2019;Zhan & Gregg, 2019) and others utilize an elasto-plastic or visco-elasto-plastic rheology (e.g., Gerbault et al., 2012Gerbault et al., , 2018Novoa et al., 2022;Souche et al., 2019). However, thermal processes, especially volume changes due to thermal expansion are rarely considered. Studies which do consider volume change due to thermal expansion are limited to a purely elastic rheology, neglecting viscous or plastic deformation of the host rock (e.g., Furuya, 2005;Kohsmann & Mitchell, 1986;Wang & Aoki, 2019).
To first order, thermal stresses can be estimated by taking the mechanical equation of state (e.g., Turcotte & Schubert, 2014) where the scalar values of pressure and density change are related to the trace of the stress and strain rate tensors P = −σ kk /3 and d ∕ = −̇d (repeated indices imply summation). Assuming an isochoric process (i.e., constant volume) and expressing dP This shows that thermal pressurization is linearly proportional to the temperature change with the ratio of the thermal expansion coefficient and the compressibility being the factor of proportionality. The ratio of the thermal expansion coefficient and the compressibility in intact rocks is typically on the order of 1 MPa K −1 . Considering that the temperature difference between rapidly injected magma bodies and their host rocks can easily reach several hundred degrees we can estimate that thermal pressure change can reach several hundred MPa. Moreover, in case of partially molten rocks, the volume change of melting/crystallization should be considered as well, implying even larger pressure changes. Based on these simple estimates, thermal pressurization can potentially exceed a near-lithostatic background pressure, potentially reaching the brittle yield stress in the host rock or significantly impacting the pressure-temperature (P-T) conditions in the magma body. Therefore, it appears feasible that thermal expansion related stresses can generate significant pressure and stress anomalies, that might even lead to thermal cracking around rapidly emplaced, cooling upper crustal magma bodies.
Our aim in this paper is to quantify stresses and deformation generated by thermal expansion/contraction around cooling magma or mush chambers in a visco-elasto-viscoplastic host rock. To do so, we have developed a new numerical code that can be used to quantify volume changes due to elastic compressibility, thermal expansion, and plastic dilation in a thermodynamically consistent manner. Besides that, the plasticity model we use considers both shear and tensile failure. Since we focus on isolating and quantifying the effects of thermal stresses around magma or mush chambers, we exclude other processes from our models. Hence we consider single phase flow (i.e., no phase separation), constant material parameters that are typical of intact granites, no background tectonic stresses. Also, we assume a pre-existing magma body (i.e., we do not model the emplacement mechanism), where the magma body has an initially elevated temperature (and thus lower viscosity and density), but otherwise is identical to the host rock. We carry out 2D plane strain thermo-mechanical simulations applied to a magma chamber with a horizontally prolate ellipsoidal geometry. We are using rheological models of increasing complexity to show the difference between a purely elastic, a visco-elastic or a visco-elasto-plastic rheology. Furthermore, we compare the influence of thermal stresses and visco-elasto-plastic relaxation (without thermal expansion) on the stress evolution around cooling magma chambers. Finally, we discuss the potential roles that KISS ET AL.

10.1029/2022JB025341
3 of 25 thermal stresses and thermal cracking might have on the evolution of magmatic plumbing systems and on the evolution of magmatic-hydrothermal systems around plutonic bodies. Our results highlight the importance of considering thermal stresses to quantify deformation and fracturing around magma chambers, when time scales over a thousand years are considered.

Governing System of Equations
We assume slow (i.e., negligible inertial forces), compressible, single velocity (i.e., multiple phases may be present, but phase separation is excluded), visco-elasto-viscoplastic deformation. The governing system of equations in 3D is where Equations 3-5 have been derived from the conservation of mass, momentum, and energy respectively. Equations 6 and 7 are constitutive relationships between volumetric and symmetric-deviatoric components of stress and strain rate tensors (Schubert et al., 2001). Indices ijk correspond to coordinate axes 1, 2, and 3 and repeated indices imply summation (Einstein notation). The strain rate tensor (̇= ) can be decomposed into a volumetric part . Our rheological model features a viscous, an elastic and a viscoplastic element in a Maxwell-type coupling for shear deformation and a thermo-elastic, and a viscoplastic element in a Maxwell-type coupling for volumetric deformation ( Figure 1, see Table 1 for parameters). This formulation accounts for processes such as compressibility, thermal expansion, plastic dilation, force balance, adiabatic heating, heat conduction, heat production due to dissipative deformation and radioactive heating, in a thermodynamically self-consistent way (for detailed derivation see Appendix A). It is worth noting that the interplay between the aforementioned processes results in a non-linear behavior which is further enhanced by the non-linear visco-elasto-viscoplastic rheology of the host rocks.

Numerical Implementation
Here, we present a 2D implementation of Equations 3-7, assuming plane strain conditions (i.e., component 3 of velocity is zero, and component 3 of all gradients are zero). The system of non-linear equations (Equation 1-5) is discretized on a regular Cartesian staggered grid using finite differences. The problem is solved by a pseudo-transient iteration or relaxation scheme (Räss et al., 2022;Versteeg & Malalasekra, 2007). Our implementation is a natural extension of the methods presented by Duretz et al. (2019) and Kiss et al. (2019) to resolve thermo-mechanical coupling for incompressible, purely viscous materials. However, we consider a non-linear visco-elasto-viscoplastic rheology, which is why we introduce new internal variables (i.e., stresses are split into trial stresses and viscoplastic stress corrections). We chose P tr , v i and T as the primary variables, and as a result Equations 4-6 are recasted in the following form: where are derivatives with respect to pseudo time ω, and are integrated in an explicit, forward Euler manner.
The terms can also be regarded as residuals of the conservation equations, decreasing during the iteration cycle. Superscript it−1 denotes values from the previous iteration and old denotes a fully converged value from the previous time step accounting for semi-Lagrangian advection. Therefore the total time derivates denote According to the small strain formulation, we neglect the rotational terms in the time derivative of the stress tensor. The last term on the right hand side of Equation 9 is introduced to dampen oscillations of the momentum residuals and hence accelerate convergence. In addition, viscosity, stress, and density are updated in an iterative manner as: To improve convergence and robustness, we employ a logarithmic relaxation scheme on the effective viscosity. In our staggered grid discretization the non-diagonal components of τ ij , tr and dev,ve are located in the vertices. Therefore the effective viscosity η ve is calculated not only at cell centers, but also at the vertices, using interpolated values of dev,vis II and T. To add the plastic correction to the non-diagonal components of τ ij , interpolated values of tr II , tr II and η ve calculated at the vertices are used.
For each physical time step Equations 8-17 are iterated until the residuals (left hand side) of Equations 8-10 reach a given tolerance value (respectively set to 10 −17 s −1 , 10 3 Pa/dy and 10 −3 K/dt in infinity norm). In addition to checking for convergence of the momentum Equation 9, we check the residuals of the additive strain rate decomposition (Equation 7) as well. This ensures that a solution of the local nonlinear problem has been found. At this point, a fully implicit solution, equivalent to backward Euler time discretization, is achieved and all non-linearities are converged.

Viscoplastic Return Mapping
The importance of viscoplastic regularization in geodynamic applications has been extensively discussed by de Borst and Duretz (2020) and Duretz et al. (2020). In essence, a viscoplastic formulation alleviates the problems associated with rate-independent plasticity (i.e., mesh dependence) and improves convergence. The implementation presented by Duretz et al. (2020) is based on the formulation of Perzyna (1966), where viscoplastic regularization is achieved by a priori fixing a viscosity value ( Figure 1, η vpl ). This kind of regularization is straightforward to implement for a linear yield function. However, we consider a piece-wise linear yield function (F) and flow potential (Q) to account for volumetric plastic strains. We have found that the equivalent formulation of Duvaut and Lions (1972) is more straightforward, when a characteristic relaxation time (instead of viscosity) is fixed a priori (Simo et al., 1988). Besides its simplicity, this implementation has the benefit of producing a uniform overstress (as a function of the distance from the yield along the return map) for all segments of a non-linear yield function.
Our plastic yield function is a piece-wise linear combination of a Drucker-Prager (F DP ), a tensile (mode-1, F M1 ) and a pressure limiter yield (F PL ), considering only the dependence on the first-(i.e., mean stress σ m = −P) and second stress invariants (τ II ). The composite yield function ( Figure 2) is formulated as According to the rate-independent non-associated plastic flow rule, where the two terms on the right hand side represent the volumetric and the deviatoric components of the plastic strain rate tensor. Viscoplastic regularization is achieved by scaling rate-independent plastic strain rates with the ratio of time increment and a relaxation time (denoted by χ) (20)

of 25
Our composite yield function exhibits two corners, one of them is at the intersection of the pressure limiter and tensile yield segments and the other one is at the intersection of the tensile and the Drucker-Prager yield segments. We use a typical non-associated flow potential (Q) with dilation for the Drucker-Prager yield and associated flow potentials for the tensile and the pressure limiter yield stress. However, considering only the potential functions corresponding to the linear segments ( Figure 2, regions I, III, and V) is insufficient, and potential functions must be created for the corner regions too (Figure 2, regions II and IV). For linear yield functions and the applied plastic potential functions, the plastic multiplier can be expressed analytically, in a closed form, as described in the following pseudo-algorithm. . The region where trial stresses violate the yield is indicated by the contoured area. This area is divided into five domains, where different plastic flow potentials are defined, corresponding to the three linear segments and the two corner regions. Return mapping in the tr − tr II plane happens orthogonally to the colored contours. However, the angle of return mapping and the domain boundaries shift as a function of the ratio of η ve and dt/β as shown for a ratio of 1 in panel (a) and for a ratio of 3 in panel (b). In this figure Q I is enlarged for better visibility, as it would be barely visible otherwise.

of 25
The Drucker-Prager and tensile yield functions and the corresponding plastic potentials are often used for geodynamic applications in an identical form (e.g., Duretz et al., 2021;Rozhko et al., 2007). The corner domains and the corresponding plastic potential are defined according to Drucker's postulate (Drucker, 1952). As we use a strain rate driven formulation, we avoid any potential issues arising from the non-unique total stress to plastic strain rate relationship in the corner domain (Ottosen & Ristinmaa, 1996).

Reference Configuration
All simulations presented here are two dimensional, plane strain, applied to a prolate ellipsoidal magma body with its long axis perpendicular to the 2D cross section. Regarding the initial temperature field, we explore two end-member cases. In the first case the magma chamber is represented as a sharp temperature anomaly and in the second case the magma/mush chamber is represented as a smooth temperature anomaly (Figures 3a and 3b, respectively). The first end-member with the sharp temperature anomaly could represent a rapidly formed magma body that did not have sufficient time to cool. On the other hand, the second end-member is representative of a long lived magmatic system. Our reference models are based on a 10.5 (sharp anomaly) and 17 wide (smooth anomaly) and 10.5 km deep model domain, with a flat initial topography and 2 km of sticky air (low density, low viscosity layer) on top. We use a free surface boundary condition on top and fixed free slip conditions on the other boundaries. The used model configuration ensures that boundary effects are small ( Figure S1 in Supporting Information S1). We apply a constant 20°C in the sticky air layer and a constant 450°C at the bottom boundary. The side boundaries are insulating (i.e., zero heat flux). The initial, background temperature field is the equilibrium geotherm, resulting from the boundary conditions, a constant thermal conductivity (3 W m −1 K −1 ) and a constant 9 of 25 radiogenic heat production rate (10 −6 W/m 3 ). The magmatic intrusion is implemented as a circular high temperature (750°C) domain, with a radius of 1.5 km and center at 5 km depth. A corresponding (in 3D) prolate ellipsoid with the semi minor axes of 1.5 km and aspect ratio of 1:4 has a volume of ca. 57 km 3 . Such magma volumes are in agreement with the estimated volumes of individual intrusions in the Torres del Paine intrusion complex (Leuthold et al., 2012). The initial stress field and the corresponding density field are calculated using a (temporally) isothermal, purely viscous Stokes solution. Buoyancy stresses in this configuration are negligible (∼0.2 MPa), hence the resulting stress field is nearly lithostatic. The input parameters and material parameters are defined as listed in Table 2, unless specifically stated otherwise.

The Purely Thermoelastic Case
As a reference, we present results from a model that considers a purely elastic rheology. We consider our reference configuration (Figure 3a) with constant parameters of α = 3 × 10 −5 K −1 , β = 10 −11 Pa −1 , and μ = 6 × 10 10 Pa (giving a Poisson's ratio of 0.25), which are typical values for intact granite. The general behavior of the system is illustrated in Figure 4. One can observe that the temperature change is largest at the contact of the magma body and its host, and it gradually decays with increasing distance from the thermal anomaly. As a result of thermal expansion/contraction, pressure changes are observed that are linearly proportional to the temperature change. However, due to the non-zero volumetric deformation, the magnitude of thermal pressurization is about half of what is expected based on the isochoric assumption. Unlike thermal pressurization that can be positive or negative, the second stress invariant is proportional to the absolute value of temperature change. For a purely thermoelastic case, the factor of proportionality is largely time independent due to the lack of stress relaxation mechanisms. Finally, total displacements in our models due to thermal expansion and contraction do not exceed a few meters at any point in time over the course of the entire simulation time of over 300 kyr. As a result detecting such processes using real time monitoring of surface deformation above a magma chamber is challenging.

Viscous Relaxation of Thermal Stresses
To illustrate the effects of viscous relaxation, we initially considered a constant viscosity and we carried out simulations with the same material properties as in the purely thermoelastic case. In this case, the constant viscosity is included in the rheological model. The results indicate that viscous relaxation has little effect initially. However, a further increase of thermal stresses leads to the gradual decrease of the magnitude of deviatoric where q x is the horizontal conductive heat flux. For both cases, the initial stress field is near lithostatic. Since the overall size of the smooth thermal anomaly is larger, we increased the model width for configuration (b) to minimize boundary effects ( Figure S1 in Supporting Information S1).
KISS ET AL.

10.1029/2022JB025341
10 of 25 stresses ( Figure 5). Consequently, thermal stresses are not sustainable indefinitely, unlike in the purely thermoelastic case. The timescale of viscous relaxation is shorter for smaller values of viscosity in agreement to the Maxwell viscoelastic timescale. For example, considering a typical lithospheric viscosity of 10 23 Pa s, one can expect relatively insignificant viscous relaxation during the first 25 kyr. However, even if one considers a viscosity of 10 21 Pa s, which is unrealistically small for most upper crustal rheologies, thermal stresses can reach several hundred MPa, which can be sustained for about a thousand years.
In the previous simulations, we considered cases with constant viscosity in the entire model domain. However, the viscosity of magmas is significantly lower than that of the host rock. To account for the possible effects of low magma viscosity, we carried out simulations where the viscosity of the host rock was kept constant but the viscosity inside the initial magma body was set to 10 20 Pa s (the viscosity of magmas is much lower, but due to numerical reasons we must limit the maximum viscosity contrast in our model). The models show that the decreased viscosity results in a rapid relaxation of deviatoric stresses within the magma body ( Figure 6). Hence the total pressure drop inside the magma body undergoes rapid spatial homogenization instead of following the pattern of total temperature change. Nevertheless, the decreased viscosity in the initial magma body has negligible effects on the stress relaxation in the host rock.

Thermal Stresses With a Realistic Visco-Elasto-Viscoplastic Upper Crustal Rheology
As discussed in the previous section, considering typical crustal or even asthenospheric viscosities, thermal stresses can reach several hundred MPa and can be sustained for thousands or tens of thousands of years. Such stress levels in a relatively shallow, upper crustal setting likely exceed the brittle yield stress. Therefore, we carried out simulations featuring a visco-elasto-viscoplastic rheology. In these simulations, we used a combination of a Drucker-Prager and a tensile yield, as explained in Section 2.3. Unlike in the previous subsection where we used constant viscosity, we here used a temperature dependent, power-law flow law of Westerly granite   (Carter & Tsenn, 1987). The results show that thermal stresses are indeed sufficient to trigger plastic failure around the upper half of the magma body, with the extent of plastic deformation being more prominent at shallower depths (Figure 7, Movie S1). Moreover, after the magma body has cooled sufficiently, plastic failure occurs within the magma body as well. This is explained by the fact that viscosity has an Arrhenius type temperature dependence and the plastic yield is pressure dependent. Therefore, in high temperature regions viscous relaxation dominates whereas plastic relaxation dominates in low temperature regions. Pressure has a secondary effect as low confining pressure promotes plastic deformation. In the regions with viscous relaxation, deviatoric stresses vanish at a characteristic time scale (Figures 5 and 7h). In the regions where plastic deformation dominates, stresses exceeding the plastic yield are relaxed back to the yield shortly after the loading ceases. As a result, stress and pressure variations that do not exceed the plastic yield can be preserved long after the equilibration of the temperature field (Figures 7g and 7h). Ultimately, the magnitude of pressure change and deviatoric stresses are limited by the plastic yield stress. In this particular case, deviatoric stresses of 200 MPa can be reached initially in the host rock, due to the initial increase of pressure and hence yield strength. Following this, the host rock cools after its initial heating phase and thermal pressurization is reversed, decreasing the confining pressure and the yield strength. As a result, the maximum stress levels gradually decrease to around 80-100 MPa. Notably, after sufficient cooling and crystallization, as a result of thermal contraction and the related decompression of the magma body, shear and tensile failure can occur inside the intrusion (Figure 8).
Despite the relatively dynamic nature of such systems in terms of pressure and stress evolution, the finite strain and total displacements are rather small, hardly observable on the macro scale. Nevertheless elastic bending of the crust near the surface can result in tensile failure albeit under small values of finite strain.

Sensitivity to the Size and Ellipticity of the Magma Body
In order to assess the sensitivity to the size and ellipticity of the magma body, we carried out nine simulations with different initial geometries of the magma bodies. We defined the geometry using the following equation: where the center of the ellipse is at x c = 0, y c = −5 km, and the semi axes are varied independently as 0.5, 1.0, and 1.5 km. The first order effects of thermal stresses are displayed for various intrusion geometries in Figure 9. The results show that the size of the fractured volume around the intrusion is proportional to the size of the intrusion, and the shape of the fractured zone is similar to the shape of the magma body.

Multiple Pulses in a Magma/Mush Chamber
Magmatic systems are generally thought to evolve incrementally, as a result of several smaller pulses of magma that are emplaced as dikes or sills in a mushy reservoir (e.g., Cashman et al., 2017;Christopher et al., 2015;Putirka, 2017). In order to test whether our general findings holds in this case as well, we carried out a simulation featuring several magmatic pulses. The magmatic pulses are introduced as instantaneous heat pulses while mass transfer and the resulting stressing of the host rock are ignored. The temperature field is set to a uniform 750°C Unlike the previous simulations where we used a sharp thermal perturbation as an initial condition, we here use a relatively smooth thermal perturbation (Figure 10a), representing a hot mushy zone around a liquid-dominated magma chamber. Then, new elliptic intrusions are added after 2, 5,000, and 10,000 years, respectively (e.g., de Saint Blanquat et al., 2011;Zhan et al., 2012). We slightly vary the location of the recurrent heat pulses as it has been suggested by geodetic observations (e.g., Delgado, 2021). All other input parameters are identical to those of the reference model. The most important result is that significant thermal stresses develop around the new intrusions only where the temperature difference compared to the surrounding volume is sufficiently large (Figures 10b and 10c, Movie S2). Although the individual pulses may cause localized small scale thermal stresses and deformation, the overall evolution is mostly controlled by the cooling and contraction of the entire thermal anomaly as a whole (Figure 10c). Because of that, the final state of volumetric plastic strain is similar to that of the reference model apart from some minor perturbations (Figure 10d). This implies that for natural magmatic systems, it is the accumulated thermal anomaly of many pulses i.e., of key importance for the development of thermal stresses.

Sensitivity to Stress History and Comparing Thermal Stresses to Stresses Induced by Melt Intrusion
In this section, we estimate the sensitivity of our results to the stress history. We use our reference configuration with the diffuse temperature anomaly (Figure 3b) and with a compressibility of 2.1 × 10 −11 Pa −1 . As a reference, we first model the evolution of thermal stresses in this model by letting it cool. A notable difference compared to the previous simulations is that due to the initially smooth temperature field the temperature evolution and hence the build up of thermal stresses is slower. Furthermore, the evolution of the system is dominated by cooling and hence contraction in this case. As a result the maximum magnitude of stresses is lower, on the order of 70 MPa. Despite the lower level of stresses, thermal cracking still takes place, although with smaller intensity in the vicinity of the magma body.
To quantify the sensitivity to stress history we compare these reference results with results of simulations including an initial pressurization of the magma chamber ( Figure 11). In the beginning we pressurize the magma chamber by modeling the gradual injection of additional magma, introducing a source term in the mass balance equation. We do not model the transport mechanism, instead we restrict ourselves in quantifying the stress evolution due to the injection in a mechanically confined volume. We stop the injection after a maximum pressure change (compared to the initial pressure) of 25, 50, and 75 MPa is reached. In our configuration, non of these injection events result in fractures that connect the magma body with the surface. As the injection stops the pressure in the middle of the magma body starts to drop until a quasi-steady state is reached at about 70 MPa below the starting value (Figure 11d). The initial pressurization has little influence on the value of this quasi-steady state pressure, apart from an increase in the time needed to reach it. In order to see the effect of thermoelastic stresses, we have repeated the same simulations, using identical parameters, but we disabled the effects of thermal expansion/contraction in the mechanical problem formulation. In these simulations, we see the effects of visco-elasto-plastic relaxation after the initial pressurization stops. However, the pressure in the center quickly reaches a steady-state value that is larger than the initial value (Figure 11d). Comparing these results with the fully coupled results shows that on short timescales (<1 kyr) visco-elasto-plastic deformation dominates the stress evolution, but thermal expansion and contraction become increasingly important at longer timescales.

Limitations Due To Simplifying Assumptions
We have presented results from numerical simulations that show the effects of thermo-elastic strains on the stress evolution of cooling magmatic bodies. In our treatment, we have neglected processes that are potentially important to quantify thermal stresses and the development of fracture networks or dikes around cooling magma bodies. On one hand, our simplifications constitute a step forward toward a more complete formulation. On the other hand, these simplifying assumptions are useful, because they allow us to isolate the effects of thermal expansion/contraction from other processes.
Most notably, we have focused on a single-phase formulation and did not include a percolating fluid or magma phase, which could reduce the effective confining pressure and hence promote localized plastic failure, by channelized porous flow (e.g., Katz, 2008;Keller et al., 2013;Schmeling et al., 2019). Another assumption we made was to neglect the volume change and latent heat of melting and crystallization, and used uniform thermodynamic parameters. By making the first assumption (i.e., single-phase flow), we decreased the likelihood of failure and the magnitude of plastic strain. In addition by considering uniform thermo-elastic properties, we underestimate the total volume change. Consequently, the results presented in this paper should be treated as a lower bound on the extent of fracturing around cooling magma bodies.

The Role of Thermal Stresses During the Evolution of Magmatic Plumbing Systems
Based on our model results, we can asses that thermal stresses likely cause a pressure change on the order of 100 potentially reaching 200 MPa. Such stress levels are comparable to the value of background pressure in the upper crust, at about 5 km depth. These pressure anomalies are accompanied by deviatoric stress anomalies of a

10.1029/2022JB025341
17 of 25 similar magnitude. The magnitude of these stress anomalies is limited by the yield strength, while their duration is controlled by viscous relaxation. For magmatic bodies that are occurring at slightly deeper levels, the yield strength is higher due to the higher confining pressure, and therefore larger pressure and stress anomalies may develop. However, due to the downward increasing temperature, the temperature difference between a magma body and its host rock is decreasing with depth. The reducing temperature difference results in a decreasing thermal pressurization after a certain depth level is exceeded. The magnitude and the distribution of thermal stresses is controlled by characteristic scales of temperature change due to heat conduction. Accordingly, the affected volume increases over time. Based on our simulations, thermal stresses start to dominate overall after about 5-10 kyr, while stress changes on shorter time scales are likely related to magma transport. On the time scale of activity of magmatic plumbing systems, thermal stresses may play a significant role.
Thermal stresses perturb the background stress field. As dykes and sills are directed by the principal stress trajectories, thermal stresses may play a significant role in the orientation and location of new intrusions (e.g., Maccaferri et al., 2011).
Besides, the direct influence thermally induced stress perturbations can result in thermal cracking and in the formation of a fracture network around the magmatic body. Some of these fractures may develop into dykes as new pulses of magma arrive from a deeper source, presenting a potential to control the evolution of the plumbing system.
Despite the relatively large values of stress perturbations, the resulting strains and displacements are rather minor compared to what can be observed by field mapping or by monitoring active surface deformation. To demonstrate this point we traced a chain of passive markers that were originally located horizontally at 3,200 m depth. The maximum displacement is −7 m directly above the center of the magmatic body. The displacement magnitudes gradually diminish as the distance to the center of the magmatic body increases. By tracing a similar marker chain Figure 11. Pressure and temperature fields (a), deviatoric stresses (b), and plastic volumetric strain (c) using our reference model setup with a diffused initial temperature anomaly (Figure 3b). Although the model domain in wider, in panels (a-c) we zoom in to the same view as in previous figures. In panel (d) the pressure evolution in the center of the magma body (0, −5 km) is compared for different magnitudes of initial pressurization. The results of no initial pressurization (ΔP inj = 0) are the same as those displayed in panels (a-c). Furthermore, results of fully coupled thermo-mechanical simulations (solid lines) including thermal expansion and contraction are compared with results of simulations without thermal expansion and contraction (dashed lines).
KISS ET AL.

10.1029/2022JB025341
18 of 25 at the surface, we see less than 5 m displacement that takes place over more than 300 kyr, resulting in negligible deformation over observational timescales.

The Manifestation of Thermal Cracking in the Field
Our model results suggest that thermal expansion and contraction have a significant effect on the stress state of a magma body and of its host rock. Stresses induced by thermal expansion and contraction are sufficient to trigger brittle failure around cooling magma bodies in a shallow, upper crustal setting. However, our models are based on continuum mechanics and a continuum theory of plasticity therefore we cannot resolve individual fractures and their characteristics cannot be directly obtained. Nevertheless, the plastic strain which results from our model can be interpreted as proxy for fracture density and also as a proxy for porosity if volumetric plastic deformation is considered. Accordingly, plastic volumetric strain is rather small (less than 0.6%) and it is not strongly localized.
Using simple estimates, that would approximately translate into a single 20 cm wide dike or into 1,000 2 mm wide joints in every 41 m of host rock (which is the grid spacing used). Consequently, plastic deformation predicted by our models may manifest on the field as a few prominent dykes or veins, or as a set of numerous fine joints, similar to exfoliation joints in granites or columnar joints in basalts, or some combinations of the two.

Thermal Cracking During the Development of Magmatic-Hydrothermal Systems
Based on the model results, we hypothesize that thermal stresses might play a considerable role during the development of magmatic-hydrothermal systems.
First, thermal cracking results in the development of fractures and joints around a cooling magma body. The volume affected by thermal cracking can extend several km away from the original magma-host contact, mostly above the magma body. This fractured volume can act as a permeable fluid pathway, which might enable or enhance the development of hydrothermal circulation around the magma body and chemical exchange between the fluids and the magma (e.g., Ruz-Ginouves et al., 2021). Moreover, as the magma body cools and crystallizes, fractures or joints may form inside the original magma volume, which can enable fluid circulation inside the crystallizing, but still relatively warm, plutonic body. The presence of such conditions might be necessary (but not sufficient) for segregation processes to take place and to leach metals from the fresh igneous rock, and thus presents a potential source of mineralization.
Second, due to thermal contraction in a relatively well confined and closed system, cooling magma bodies undergo decompression even when the magma body remains essentially stationary ( Figure 12). This is of potential importance as the solubility of volatiles in melts is primarily a function of pressure. To illustrate this, we used the water solubility models from Volatilecalc, presented by Newman and Lowenstern (2002) for rhyolitic melts ( Figure 12). Therefore, if such a plutonic body has sufficient amounts of volatiles and it undergoes decompression due to cooling, volatiles might be expelled from the melt and appear as a free phase. Thus, thermal contraction induced decompression might introduce an additional fluid source for the magmatic-hydrothermal system.

Conclusions
We presented a numerical method that is suitable to quantify stress evolution related to thermal expansion/ contraction in an upper-crustal setting with visco-elasto-viscoplastic rheologies including both shear and tensile failure.
Our results demonstrate that thermal stresses around upper crustal magma bodies are significant as stress anomalies can reach or even exceed the background lithostatic pressure. Pressure anomalies are proportional to the temperature change, but viscous or plastic relaxation might limit their magnitude or duration. The host rock nearby the magma body experiences significant pressurization upon heating. At the same time, cooling and thermal contraction causes significant decompression in the magma body.
Moreover, thermal stresses are likely sufficient to create an extensive fracture network around an upper crustal intrusion by brittle failure. The exact depth where brittle failure may occur is dependent on the rheology of the rock and on the depth of the magma body.
Over the scale of several to 100 kyr, thermal stresses might contribute to the development of the magmatic plumbing system as pressure perturbations and the developing fracture network might influence the location 19 of 25 of new intrusions. Furthermore, we speculate that the appearance of a volatile phase and the development of a fracture network around the magmatic bodies has the potential to one of the main processes that control magmatic-hydrothermal alteration around magmatic bodies. Hence, thermal stresses may play an important role during ore mineralization or post-volcanic activity as well.

Appendix A: Thermodynamic Admissibility of the Governing Equations
In this appendix we show the thermodynamic admissibility and consistency of the governing equations, based on classical irreversible thermodynamics (e.g., De Groot & Mazur, 2013;Müller & Müller, 2009), which has previously also been used in a geodynamic context to demonstrate the thermodynamic admissibility of two phase flow formulations (Yarushina & Podladchikov, 2015). For the description of recoverable and dissipative bulk and shear deformation we follow Landau and Lifshitz (2013). We assume single velocity deformation in isotropic, visco-elasto-plastic, compressible materials with constant chemical composition.

A1. Local Thermodynamic Equilibrium
We use the local thermodynamic equilibrium (LTE) assumption to relate equilibrium thermodynamic relationships to continuum mechanics. In essence, the LTE states that equilibrium thermodynamic relationships are applicable locally and instantaneously, even if the system is not in global equilibrium (e.g., there are pressure or temperature gradients). This can be utilized to formulate a relationship between the increment of specific total energy (dE) and the sum of the increments of heat (TdS), kinetic energy, potential energy, elastic strain energy, and nuclear energy (radiogenic heating): Here, elastic strain rate denotes all recoverable deformation, including thermal expansion or contraction. We consider an isotropic, Maxwell visco-elasto-viscoplastic rheology, which implies an additive strain rate decomposition and uniform stress on each rheological element. In the limit of purely hydrostatic conditions and entirely recoverable volumetric deformation, the elastic strain energy equals −Pdρ −1 (e.g., Landau & Lifshitz, 2013;Müller & Müller, 2009).

A2. Local Balance Equations
The local balance equations of mass, linear momentum, energy, and entropy in a Lagrangian form are given as follows: where ρ −1 , E and S are respectively specific volume, specific total energy, and specific entropy. , and , are respectively the non-advective specific momentum, specific energy, and specific entropy fluxes, defined as Thermodynamic admissibility is ensured if the source of specific entropy, Q S , is non-negative, which we will demonstrate in the following sections.

A3. Solving for Q S
To relate this LTE to the balance equations, we express the LTE (Equation A1) using increments with respect to time and multiply it by ρ: Note that we have applied the chain rule to simplify the kinetic energy term (second term on the right hand side). Then we substitute Equations A6-A8 into Equation A9 to replace the time derivatives and solve for TQ S After substituting the momentum flux (Equation A6) into the third term on the right hand side and using the sum rule, the potential energy cancels out Now substituting the energy flux (Equation A7) into the first term on the right hand side, using the difference and the product rules, the following terms remain KISS ET AL.
Finally, substituting the entropy flux (Equation A8) in the right hand side and simplifying yields

A4. Demonstrating the Non-Negativity of Q S
Non-negativity of entropy production (Equation A14) is guaranteed if all terms on the right hand side are non-negative (T [K] is non-negative). It is easy to see that dissipation due to heat conduction and radioactive heating, the first and the third term on the right hand side, respectively, are guaranteed to be non-negative for any non-negative λ and Q r . Showing the non-negativity of the second term on the right hand side (dissipative work), however, requires to explore the rheological models.
The strain rate (or velocity gradient) tensor can be expressed as a sum of its volumetric, symmetric-deviatoric, and antisymmetric parts while the stress tensor is symmetric (e.g., Landau & Lifshitz, 2013) and can be decomposed into a volumetric (hydrostatic) and a deviatoric part where = − 3 is thermodynamic pressure. Our rheological model for shear deformation is based on the Maxwell (serial) coupling of a viscous an elastic and a viscoplastic element. The rheological model for volumetric deformation consist of an elastic element, that represents the pressure-volume-temperature equation of state in an incremental form, in a Maxwell coupling with a viscoplastic element. The relationships between stresses and strain rates, broken down for deviatoric-symmetric and volumetric parts, are:̇d ev = 2 ⏟⏟⏟ dev,vis Considering that for any two second order tensors, therefore (̇−̇e l ) =̇d ev,vis +̇d ev,pl −̇v ol,pl .
After substituting the stress-strain rate relationships from Equations A17, A18, and A20 becomes For any positive η, viscous dissipation is non-negative, and plastic dissipation is zero if the yield is not reached (since ̇= 0 ). Moreover, if a rate-independent plasticity model is admissible (χ = 1), the same plasticity model with Duvaut-Lions regularization must be admissible too (since 0 < χ ≤ 1). To show the general admissibility, we demonstrate the admissibility of each of the five plasticity models in the rate-independent limit, that were introduced in Section 2.3.
Plastic dissipation for the pressure limiter yield (domain I) after substituting derivatives of Q I and expressing P from F PL = 0 iṡp which is non-negative as long as (σ T − δσ T ≥ 0) since ̇ is never negative.
Plastic dissipation for the tensile yield (domain III) after substituting derivatives of Q III and expressing P from which is non-negative as long as (σ T ≥ 0).
Plastic dissipation for the Drucker-Prager yield (domain V) after substituting derivatives of Q V iṡ pl * =̇( II − sin ).
It can be shown that τ II − P sin ψ ≥ 0 by taking F DP = τ II − P sin φ − C cos φ = 0. Since τ II and C cos φ are both positive, τ II − P sin ψ ≥ 0 as long as 0 ≤ ψ ≤ φ, because the sinus function strictly monotonically increases from 0° to 90°.
The plasticity model for the corner regions (II and IV) can be generalized, as the only difference is the pressure-stress coordinates of the corners P = P C and τ II = τ C . Plastic dissipation after substituting derivatives of Q II and substituting P = P C and τ II = τ C iṡ which is guaranteed to be non-negative if the numerator is non-negative. The numerator can be reformulated as Any trial stresses in the corner regions can be expressed in the following form, where 0 ≤ ϒ ≤ 1 for the first corner and 1 ≤ Υ ≤ 1 sin for the second. Substituting Equation A27 into the numerator Equation A26 and simplifying results in which is guaranteed to be non-negative for P C ≤ 0 (first corner, domain II) and for any τ C ≥ P C if ϒ ≥ 1 (second corner, domain IV).

A5. The System of Governing Equations in Their Final Form
The governing equations used in the manuscript (Equations 3-7) are all based on the balance laws (Equations A2-A5), fluxes (Equations A6-A8) and constitutive relations (Equations A17-A18).
The divergence of the flux and the conductive dissipation (first and second terms on the right hand side) can be merged, using the product rule, The same equation can be obtained by substituting the LTE (Equation A9) and the fluxes (Equation A7) into the conservation of energy (Equation A4).
So far, we have three equations and six unknowns (ρ, v i , S, P, τ ij , T). In order to close the system of equations, we first formulate entropy increments as a function P and T d = The definitions of C P , α, and a Maxwell relationship can be used to express the coefficients of dT and dP gives which is identical to Equation 5. By expressing entropy as function of pressure and temperature we thus reduced the number of unknowns to 5 and by including the constitutive relationships for bulk and shear rheology (Equations 6 and 7) we obtain a closed system of equations, that is both thermodynamically admissible and self-consistent.

Data Availability Statement
We have developed a julia code to solve the governing equations. The full source code to reproduce the reference simulation (Figure 7) is available under https://zenodo.org/record/6958273 (https://doi.org/10.5281/ zenodo.6958273). The other simulations can be reproduced by modifying the reference case with the parameters described in the manuscript.