Thermo‐Hydro‐Mechanical Model and Caprock Deformation Explain the Onset of an Ongoing Seismo‐Volcanic Unrest

Modeling seismicity at volcanoes remains challenging as the processes that control seismic energy release due to fluid transport, heat flow, and rock deformation are firmly coupled in complex geological media. Here, we couple fluid‐flow and mechanical (deformation) simulators (TOUGHREACT–FLAC3D) to reproduce fluid‐induced seismicity at Campi Flegrei caldera (southern Italy) in isothermal (HM) and nonisothermal (THM) conditions. The unique ability of the Campi Flegrei caprock to withstand stress induced by hot‐water injections is included in the model parametrization. After pore pressure accumulation is guided by a combination of thermal and hydromechanical interactions, contrasting compressive and extensional forces act on the basal and top parts of the caprock, respectively. Then, pressure perturbation and caprock deformation induce fractures that allow hot fluids uprising to pressurize the overlying fault, driving it toward failure and triggering seismicity. Under similar mechanical boundary conditions, the induced thermal effects prompt seismic slip earlier but with higher seismic magnitudes when (1) thermal equilibrium is preserved and (2) the thermal contrast is enhanced due to increased fluid injection temperatures. The results indicate that numerical models of volcano seismicity must consider the influence of rock‐sealing formations to obtain more robust, accurate, and realistic seismic predictions at volcanoes. The proposed models satisfactorily reproduce the magnitude–depth distribution of the swarm (October 5, 2019), preceding the two strongest earthquakes recorded in 35 years at the caldera (3.1 and 3.3—on December 6, 2019, and April 26, 2020, respectively) using hot‐water injection from depth.

, but little attention has been so far paid to the reproducibility of the observed seismicity quantitatively using numerical simulation approaches. Yet, numerical modeling of seismicity could enhance volcano monitoring and forecasting, improving volcanic risk assessment and hazard mitigation.
This level of detail on volcanic structures makes CFc the ideal laboratory to simulate the seismic response to fluid injections during volcanic unrests. Between September 1983 and April 1984, it was the increase in number and intensity (up to 3.8 Mw) of shallow (<4.5 km depth) earthquakes in the caldera that convinced the authorities to evacuate >40,000 people from the city of Pozzuoli (Kilburn et al., 2017). The increase in seismicity was not followed by an eruption. There is wide agreement that this unrest was induced by a relatively shallow magma reservoir producing magmatic intrusions, imaged between 3.5 and 4.5 km by seismic attenuation tomography in the same area where it is modeled by deformation studies (Amoruso, Crescentini The presence of shallow magma activity (depths of 2-4.5 km) feeding the hydrothermal systems during recent unrests (from 1985 to 2014) is more debated by recent studies. Some still attribute the release of CO 2 -rich fluids to a shallow magmatic sill (Chiodini, Paonita, et al., 2016) while other interpret the unrest as due to the drying of the base of the hot hydrothermal system (Moretti et al., 2018-see Troise et al., 2019, for a recent review), with no melt between the surface and a deeper magmatic reservoir (∼7 km depth- Zollo et al., 2008). The repeated unrests at the caldera promote a long-term accumulation of stress in the crust, with the caldera possibly evolving toward conditions more favorable to eruption (Kilburn et al., 2017). Moreover, after almost 40 years of earthquake magnitudes never exceeding 3M, seismic energy release is finally increasing: at the time of submitting this paper, the most significant event (3.3 M) has been recorded at 2.57 km depth on April 26, 2020, under the Solfatara/Pisciarelli craters (Bollettino INGV-OV, April 26, 2020). Independently of their origin, all studies concur that hot fluid injections from depths is the dominant mechanism for the onset of microseismicity at the caldera (Chiodini, Paonita, et al., 2016;De Siena, Chiodini, et al., 2017;Moretti et al., 2018;Petrillo et al., 2013).
In this study, we carried out single-phase isothermal (HM) and nonisothermal (THM) hot-water injections by coupling multiphase multicomponent fluid and heat flow simulator (Transport Of Unsaturated Groundwater and Heat, TOUGHREACT) with a mechanical (deformation) simulator (Fast Lagrangian Analysis of Continua in three dimensions, FLAC 3D - Gan & Elsworth, 2016;Taron & Elsworth, 2009). We thus explored the THM effects produced by hot-water injections to retrieve seismic slip, and its relation to event magnitude and earthquake depth during seismic unrest. This study is thus complementary to existing studies, which have concentrated on the uplift (ground deformation or vertical displacement) magnitude and geochemical variations (Chiodini, Caliro, De Martino, et al., 2012;Chiodini, Paonita, et al., 2016;A. Coco et al., 2016;Kilburn et al., 2017). After introducing the background of our modeling approach, the modeling results are presented and discussed within the context of thermal-driven processes. Such processes are considered the main cause of the 2005-2014 unrest at the caldera based on deformation and geochemical anomalies (Chiodini, Vandemeulebrouck, et al., 2015) and therefore used as inputs for TOUGHREACT simulations of magmatic degassing (Chiodini, Paonita, et al., 2016).
The coupling between TOUGHREACT and FLAC 3D allows us to model seismic parameters, such as earthquake depths and magnitudes, in combination with deformation and poroelastic effects, an important addition considering the renewed seismic activity at the caldera. Central in this application is the inclusion of realistic heterogeneities in the subsurface discretization, especially the mechanical strength of the caprock and its role in modulating stress (Vanorio & Kanitpanyachaeron, 2015). This low-attenuation and low Vp/Vs lithological change appeared to hinder uprising magmatic fluids in 1983-1984(Akande, De Siena, et al., 2019Calò & Tramelli, 2018;De Siena, Chiodini, et al., 2017) and is still likely a major controller of unrest today. We specifically targeted the recent unrest as there is agreement that the 1982-1984 unrest was induced by shallow magmatic intrusions Amoruso et al., 2008;De Siena, Chiodini, et al., 2017;Moretti et al., 2018;Troise et al., 2019), and seismicity was very shallow and weak between these two unrests (Bianco et al., 2004), leading to inferences of changes in the rheological structure of the volcano Di Luccio, Pino, et al., 2015). The key questions for our study are (1) if the highest recorded seismicity since April 1984, in 2019 and 2020, can been induced by hot-water injections without the need of a magmatic source at shallow depths and (2) if it is better modeled by isothermal or nonisothermal conditions. Our model was constructed using structured block grids with a gradational size and local refinements that characterize the fault zone within the first 2 km depths (A. Coco et al., 2016;Troise et al., 2019).
The dimensions of the spatial model were 10 × 1 × 3 km. The parametrization of the model, discretized into six horizontal geological layers, was based on hydraulic unit information from an AGIP report (AGIP, 1987), physical characterization of the rocks (Vanorio & Kanitpanyachaeron, 2015), and tomographic studies (Akande, De Siena, et al., 2019;Calò & Tramelli, 2018;De Siena, Chiodini, et al., 2017). Grid resolutions are similar for both simulators. The maximum grid spacing along the x-and y-axes (Δx and Δy) is 500 m, while it is 100 m along the z-axis (z), giving an aspect ratio within 1:5 which ensured model accuracy and solution robustness (Abbasi et al., 2013). To reduce the computational time for simulations, the extent of the intermediate y-axis was limited to 1 km, assuming that the properties of interest are almost uniform (symmetric) along this axis. However, this simplification was necessary because it did not have any significant effects on our results. The basal aquifer is defined at 2.5-3.0 km and overlain by a caprock that is about 0.5 km thick (Calò & Tramelli, 2018) at a depth of 2.0-2.5 km. The middle and upper reservoirs are located at depths of 1.8-2.0 and 0.5-1.0 km (Akande, De Siena, et al., 2019;Petrillo et al., 2013), respectively ( Figure 1).
The explicit fault architecture includes a plane of structural weakness representing the fault with a dip of 60° and strike of 290°. The fault was modeled as 10-m thick, which is the maximum seismically resolved fault width of Bruno et al. (2017). This major fault crossed the upper discretized layers and terminated on the top of the caprock. Outcrops and available subsurface data revealed that fracture systems are generally of normal kinematics with steep attitudes, as most faults having a dip angle greater than 60°, typical of subvertical normal faults (  60° dip angle gave more robust results in this study compared to other tested dip angles (70° and 80°-Figure S2, from the supporting information), a common feature of simulation studies (e.g., Steffen et al., 2014). The fault zones were modeled as a single homogeneous entity without being partitioned into a fault core and damage zone, as the variations of properties across these zones were not the major focus in this study. The fault was represented as solid elements with ubiquitous joints oriented as weak planes, implementing the recommendations of Cappa and Rutqvist (2011) and Elsworth (2014a, 2014b). Uniform fracture spacings of 10.0 and 1.0 m were chosen for the rock matrix and fault zone, respectively, guided by the ranges reported in previous studies at the caldera-wide scale (Isaia, Vitale, Di Giuseppe, et al., 2015;Vitale & Isaia, 2014). A vertical injection source was designed and placed at a depth of 3 km at the center of the model, representing the injection of hot water from a deeper feeding system, for example, the 7-km-deep sill proposed by Zollo et al. (2008).
The initial principal stresses in the three directions (x, y, and z) and their gradients (Table 1) were determined based on tectonic reconstructions in the area (Bianco et al., 2004;Ventura & Vilardo, 2005). A stability test as well as sensitivity analyses were carried out to ensure the initial mechanical equilibrium conditions. Stress-displacement and pressure boundary conditions were applied on the model: lateral and bottom boundaries of the model are set as fixed perpendicular to the boundaries, and the Earth's surface sets as a free surface of 63.9°C mean temperature (Vilardo et al., 2015) and 0.1 MPa of atmospheric pressure. The pressure at the bottom of the model is set to 30 MPa based on a pressure gradient of 10 MPa/ km (Mandl, 2000). The injection source is defined as Neumann boundary condition with a constant mass rate, by assigning moderate "dummy elements" to the adjacent elements at the bottom of the model. The thermodynamic behavior of fluid transport with thermal energy transfer is simulated by the built-in equation of state (EOS1) module, which involves single-phase and isothermal/nonisothermal water transport. The chemical reactive transport process is not considered in this study, as the major scope aims to explore the influence of heat transfer and plastic deformation, which usually occurs in relative early intermediate    (Niemeijer et al., 2010;Yasuhara et al., 2003).
The fluid transport and concurrent heat transfer are allowed to occur in dual porosity-permeability media, which are created through the MINC-partitioning (Multiple Interacting Continua) of the primary grid with fracture spacing and corresponding aperture constrained (Barenblatt et al., 1960). Furthermore, the incorporated constitutive models allow dynamic evolution of mechanical properties, including the porosity and permeability for both matrix and fracture. A stress-dependent constitutive model is employed to quantify the fracture aperture evolution (b, in m) for the fractured medium (Min et al., 2009;Rutqvist, Wu, et al., 2002). The model is represented by an empirical relation, which is a function of nonlinear fracture stiffness (α) and applied effective normal stress (σ', in Pa) as: where σ' is the total normal stress minus pore pressure (Pa),   0 is the effective stress at which zero deformation occurs (typically, 0 Pa), b denotes the hydraulic aperture due to mechanical stress effect, b max is the maximum aperture at zero stress level (m), and b r is the residual aperture (m). The value of α was taken as 0.218 MPa -1 (Gan & Elsworth, 2014a).
We then obtained the permeability (k) evolution via Warren-Root model, a cubic flow relation as (Warren & Root, 1963): where s is the fracture spacing (m). The injection-induced volumetric strain feedback from FLAC 3D are employed to update porosity accordingly (Chin et al., 2000): where  Δ s v is the change of matrix volumetric strain,  1 m n and  n m are the matrix porosity at time t n+1 and t n , respectively.
Progressive fluid injection and accompanied overpressure could ultimately lead to modification of fault permeability owing to irreversible induced strain (Gan & Lei, 2020, and references therein), upon breaching of the modeled caprock in this study. Thus, fault shear dilation due to fault slip enhances fault permeability and porosity. The role of permeability evolution on modifying the fault effective stress regime becomes more significant when modeling thermal effect of fluid injections in geothermal systems as well as volcanoes.

Definition of Materials and Constitutive Properties
The model domain was populated with material properties, such as density, elastic moduli (bulk modulus and shear modulus), Poisson's ratio, and hydraulic properties (permeability, porosity, and fracture spacing) for the rock mass (matrix) and faulted blocks (Table 1). Mechanical properties for the fault zone as well as for the matrix were also specified, including cohesion, tensile strength, friction angle, and dilation angle ( Table 2). The fault was modeled as structural elements with zero cohesion and a lower internal friction angle than the matrix. It was also assigned a lower rigidity (shear modulus) value than the surrounding matrix, and this is reasonable in such a poroelastic medium under undrained conditions (Coco & Rice, 2002 Table 1).

Table 2
Mechanical and Thermal Properties Used for the Simulation matrix was assigned the Mohr-Coulomb model, while the fault was modeled with the strain-hardening/ softening ubiquitous joint model type of FLAC 3D .

Scenario Analyses and Simulations
The simulation conditions used in this study, such as hot-water injection rates, injection temperatures, and fault dip angle, are summarized in Table 3. Preliminary studies (e.g., Figures S2 and S3, from the supplementary appendix) revealed that a dip angle of 60° and an injection rate of 150 kg/s produced optimum results (yielding fault slips at reduced timings) in the isothermal mode, and this allowed the evaluation of the hydromechanical effect of hot-water injection. These conditions, as well as an injection temperature of 350°C and fault and caprock permeabilities of 1.0 × 10 −14 and 1.0 × 10 −18 m 2 , respectively, were used as the base case for all the simulation Scenarios A-D (Table 3).
In simulation Scenario A, the hydraulic effect of the hot-water injection rates on the timing of fault slip and seismic moment magnitudes were examined under isothermal condition. Previous research at CFc demonstrates the importance of performing nonisothermal simulations to evaluate the thermal effects of fluid injections on geochemical and deformation signals (Chiodini, Paonita, et al., 2016;Chiodini, Vandemeulebrouck, et al., 2015;A. Coco et al., 2016). While the temperature of the geological layers was maintained at a constant value, the injection temperatures (200°C, 250°C, 300°C, and 350°C) were thus varied in simulation Scenario B. At CFc, researchers have demonstrated the importance of considering permeability contrasts between geological layers (k m ) and fault zone (k f ) (e.g., Jasim et al., 2015).
By varying them as shown in Table 3, the influence of fault permeability (1.0 × 10 −13 to 1.0 × 10 −17 m 2 ) on the timing of fault slip and magnitude of seismicity was investigated in Scenario C. As noted earlier in Ta Note: Significance of bold indicates highlight the parameters (or conditions) being varied under different simulation scenarios.

Table 3 Simulation Scenarios and Conditions Used in the Study
However, we also examined the influence of caprock permeabilities in the range of 1.0 × 10 −16 to 1.0 × 10 −20 m 2 (Table 3, Scenario D).
The instability of the fault was evaluated by quantifying the evolution of Coulomb stress ratio. The shear stress (τ), normal stress (σ n ), effective normal stress (  n ), and Coulomb stress ratio (η) were obtained using the following equations (Cappa & Rutqvist, 2011;Gan & Elsworth, 2014a): In these equations, σ 1 is the maximum principal stress, σ 3 is the minimum principal stress, τ xz is the shear stress component along the x-z plane, θ is the direction of the maximum principal stress from the fault plane (Figure 2a), and P is the pore fluid pressure, which reduces the total normal stresses to their respective effective normal stresses ( Figure 2b). Thus, the expression in Equation 5 is often referred to as the effective stress law (Terzaghi, 1923).
When combining the effective stress law with the Coulomb failure criterion (Jaeger & Cook, 1979), the conditions for fault slip to occur at a fault with a specified orientation can be written as (Cappa & Rutqvist, 2011): where τ is the threshold shear stress for fault slip, c is the material cohesion strength, and μ s is the static coefficient of friction whose magnitude is defined as: where φ is the internal friction angle.
According to Baisch et al. (2010), fault slip occurs when the ratio of shear stress to effective normal stress (Equation 7) exceeds the frictional strength of the fault plane (i.e., η > μ). This ratio is also referred to as "slip tendency" (Morris et al., 1996) or "ambient stress ratio" (Cappa & Rutqvist, 2011), and it is synonymous to the Coulomb stress ratio of Gan and Elsworth (2014a). The friction angle used in this study was 30°, corresponding to a coefficient of internal friction of 0.577 (Bianco et al., 2004). Thus, seismic slip was expected when the Coulomb stress ratio was approximately 0.58. According to Collettini and Trippetta (2007), the potential of a fault to undergo slip depends on its orientation with respect to the stress field, and faults optimally oriented for reactivation are those whose shear stresses (τ given by Equation 8) greatly exceed the fault's frictional resistance. Seismic fault slips in natural earthquakes are often discontinuous owing to a slip-weakening process (Kanamori & Brodsky, 2001); therefore, we developed a friction table (   used it to model fault frictional slip behavior, and interseismic phases of frictional healing in line with the friction laws of natural earthquakes (e.g., Dieterich, 1979;Rice & Ruina, 1983;Ruina, 1980). Figure 3 is a typical member of curves previously described by Di Toro et al. (2011) for well-consolidated and noncohesive rocks, which depict exponential decay of friction coefficient with increasing slip. The authors have demonstrated that such a behavior is lithology dependent. The geology of CFc, which is characterized by pyrolcastics, tuffs and tuffites, and marine deposits especially at shallow seismogenic depths (e.g., Isaia, Vitale, Marturano, et al., 2019), could be potentially altered by more or less coexisting mechanically and chemically activated weakening mechanisms resulting in frictional strength reduction with slip enhancement (Di Toro et al., 2011).

Results
It takes a period equal or longer than 5 × 10 5 s (6 days) for building up pore pressure from the inception of hot fluid (water), before the initiation of fault slips in isothermal (hydromechanical, HM) conditions (Figure 4). In Scenario A, it was observed that for higher injection rates (1) the timing for the fault slip greatly reduces ( Figure 5a) and (2)   initial part of the plot before the fault slip peak illustrates the period of shear stress accumulation accompanied pore pressure accumulation, (b) fault mean permeability evolution, which shows the permeability enhancement that accompanied fault slip, (c) fault pore pressure evolution, showing that "pore pressure surge" coincides with the timing of fault slip, and (d) fault shear stress evolution, illustrating initial slow stress accumulation followed by rapid stress jump just before fault slip, and this accompanied the pore pressure surge.
the caprock, leading to the formation of microfractures that enhance the hydraulic properties of the system ( Figure 6) and ultimately provides conductive migration pathways for hot fluids to pressurize the overlying fault at a shallower depth. More generally, pressure builds up for a longer period (∼7.5 days) at very low caprock permeability, with a consequent delay in shear fracturing (e.g., Figures 6d, 6e, 6i, and 6j) and timing of fault slip (Figure 7). In the condition of lower caprock permeability, hydraulic and thermal pressurization are insufficient to penetrate the caprock quickly due to the dominance of lateral hot-water diffusion over vertical stressing. With an increase of permeability (1.0 × 10 −16 to 1.0 × 10 −18 m 2 ), local thermal effect combined with hydraulic pressurization to enhance the potential of shear failure (Figures 7a-7c).
In the isothermal case with caprock permeability of 1.0 × 10 −18 m 2 , the seismic fault slip at the fault surface occurs ∼7 days after the inception of injection. In the nonisothermal mode, a similar fault slip event occurs after ∼6 days ( Supplementary Information, Figure S4). These events are a consequence of the frictional instability of the fault zone when the ratio of shear stress to effective normal stress (τ/ 

Hot-Water Injections and Caprock Dynamics
Several studies have suggested a connection between fluid injections and induced microseismicity in volcanoes and have assessed it in geothermal fields (e.g., Cornet   Induced seismicity can be triggered by a thermoelastic effect, a hydromechanical (poroelastic) effect, or a combination of both (Gan & Elsworth, 2014a, 2014bKwiatek et al., 2015). A significant thermal contrast between the injected fluid and the injection formation can generate seismicity that migrates following the thermal gradient in the reservoir (Kwiatek et al., 2015). At CFc, fluids were likely injected from depth at a rate of 150 kg/s and at a temperature of 350°C, a temperature which was higher than the temperature of the injection formation (reservoir-340°C), at least until 2015 (Chiodini, Paonita, et al., 2016). The temperature difference caused a sudden expansion of the reservoir volume and prompted overpressurization in the basal reservoir (Figure 6). An average temperature of 340°C was used here to assume subcritical state for the reservoir and keep the model within the simulator's capability, also because the critical temperature was reached by an AGIP well in part of the caldera between 2.5 and 3.0 km only during the 1982-1984 unrest (AGIP, 1987;Troiano et al., 2011).
Following injection-induced effects of hot-water injection, the thermal front profile spread upon reaching the caprock, and the developed thermal stress built up acting on the caprock (Figure 6). This could lead to fracturing of the caprock due to thermal stress and/or weakening (Goodarzi et al., 2012;Preisig & Prévost, 2011). Fluids were vertically injected in this study, as is usually the case in volcanoes (Chiodini, Paonita, et al., 2016), so the vertical stress was greater than the horizontal stress via possible lateral fluid diffusion (Vilarrasa et al., 2011). This condition aids plastic strain propagation to localize, especially at the center of AKANDE ET AL.  the caprock (Figures 6 and 9). Relatively high permeability of the basal reservoir conferred on it a higher diffusivity, which allowed faster diffusion of fluid ( Figure 6). When the injection-induced pore pressure caused high-pressure build-up (overpressure) in the reservoir, it expanded the reservoir and the overlying caprock, which behaves like a plate and bends (uplifts) in response to the fluid overpressure of the hot-water injection underneath (Figure 9b, 9c, 9e, and 9f). Perturbed stresses are known to be concentrated at the caprock-reservoir interface (Orlic et al., 2011). If intensified, this condition can make the caprock susceptible to shear and/or tensile failure. The pore pressure front put the lower part of the caprock in horizontal compression, while the upper part of the caprock underwent horizontal extension which eventually yielded shear fractures (Figures 6 and 9). Concentrations of deformation at reservoir/caprock interface as observed in this study have been previously described by Vilarrasa et al. (2011) at a CO 2 -sequestration site in deep saline aquifers. As far as the caprock exists at the caldera, this deformation pattern might be expected during seismic unrests.
Besides stress state perturbation in the injection formation due to direct fluid injection (Table 3), the permeability of the caprock plays a significant role in determining the timing of seismic slip, and whether the caprock's integrity will be compromised to produce conductive pathways (Gor et al., 2013). The permeability scenarios, investigated for the caprock in relation to the seismic slip time, revealed that a caprock with a AKANDE ET AL.

10.1029/2020JB020449
14 of 26 Figure 9. Distributions of the permeability, pore pressure and shear stress in the caprock near fault slips (∼6 days after the injection started) illustrated in Figure 6, panels (a-c) Isothermal, panels (d-f) Non-isothermal, and panel (g) is the 2-D model showing the line of section A-A′ along which panels (a-f) were drawn. There is an increase in the caprock permeability by 5-10 orders of magnitude in both simulation conditions (panels a and d), depending on the initial permeability. lower permeability (e.g., 10 −19 and 10 −20 m 2 ) requires a longer time (reaching ∼8 days) for pore pressure to reach seismic slip. On the contrary, relatively higher caprock permeabilities (e.g., 10 −17 and 10 −18 m 2 ) result in seismic slip occurring sooner (Figures 6 and 7).

Fault Slips: Roles of Fault Permeability and Its Contrast With Matrix Permeability
While the accumulation of shear stress can take from months to year(s) at crustal scale before fault rupture, the results of this study show that shear stress accumulation due to fluid injections takes days (∼6 days) at CFc before seismic slip, resulting in seismic events (Figures 4 and 5, and Supporting Information, Figure S4). Such stress accumulations are often observed just before volcanic eruptions and earthquakes to the point where fault criticality is reached and the shear stress cannot be sustained any longer, leading to rock fracturing (Cappa & Rutqvist, 2011;Crampin et al., 2003;Volti et al., 2003aVolti et al., , 2003b. This resulting seismic slip characterizes both the past (Aster & Meyer, 1988;De Siena, Amoruso, et al., 2017;De Siena, Chiodini, et al., 2017)  . The seismic slips occur much earlier at higher hot water injection temperature cases as thermal effects create frictional weakening, and the slip magnitude is also enhanced at higher temperature contrasts (see the text for details). this period of hot-water injection. The shear fracturing and shear dilation were accompanied by a concomitant increase in mean permeability (10 −16 < k < 10 −10 m 2 ), which decays toward the fault's top (Figure 11). Figure 11 reveals that the mean permeability along the fault fell in the range of "seismogenic permeability" (10 −16 < k < 10 −14 m 2 ) from migrating microseismicity commonly induced by fluid injection at relatively shallow crustal depths (Talwani & Acree, 1984).
Pore pressurization of the fault zone is synonymous to reducing its shear strength (Mazzoldi et al., 2012). Thus, our simulation results in Scenario C indicate that the timing of fault slip is independent of fault permeability over the range of values considered (k f = 10 −13 to 10 −17 m 2 ). However, the magnitude of the slip distance is found to be higher at lower fault permeabilities (e.g., k f = 1.0 × 10 −16 and 1.0 × 10 −17 m 2 ), and this increment of slip distance may be due to rapid pressure build-up due to lack of fluid diffusion to the surroundings. On the other hand, exchange of fluid with the surrounding matrix prevented high-permeability faults from garnering enough pressure to produce much displacements during rupture. These results contrast with the findings of Mazzoldi et al. (2012) who reported, in a similar study, that there was no significant impact of permeability on fault slip, but variations in fault permeability determine the timing of fault rupture, which we also intuitively expected in our results. The timing and magnitude of fault slip in this study may be explained by a two-step diffusion process. The first diffusion (fluid migration) between the underlying basal reservoir (injection formation) and the impervious caprock determined the failure of AKANDE ET AL. sealings by caprocks, and the second diffusion between the fault and neighboring matrix determined the timing and magnitude of fault activation.

Comparison of HM and THM Effects
Fault reactivation generally occurs earlier under nonisothermal (THM) conditions, compared to the decoupled hydromechanical isothermal scenario (compare Figures 6f-6j with Figures 6a-6e). As the temperature difference increased in the nonisothermal mode, the time required for corresponding fault slip reduced. This implies that the thermal effect may further influence the patterns of seismic slip under the same mechanical loading condition and material properties as in the isothermal mode ( Figure 6 and Supporting Information, Figure S4). The enhanced pore pressure, reducing the effective normal stress to the point of failure criticality, was a consequence of the rate of temperature changes between the injecting fluid and injection formation (Pinyol et al., 2018) and feedback from heat energy transported by the injected fluid. The pore pressure and temperature modified the stress paths and brought the rock formation to its failure point faster.
Based on our observations, the seismicity at CFc is induced by both hydromechanical and thermoelastic effects. The thermal front opened up fractures faster, augmenting the pore pressure front which subjects the upper part of the caprock to extensional forces and fractures (Figures 6 and 9). Similarly to this work, Lu (2018) found that thermoelastic stress and fluid pressure changes acted upon partially open or hydrothermally altered fractures could enhance rock permeability during nonisothermal fluid injection. Bonafede (1991), in his seminal thermal-driven advection work, highlighted the necessary conditions for a connection between two pervious reservoirs (as the case of basal injection formation and overlying fault in this study) in a volcanic setting. Observations have shown that the heated deep fluid rising from depth might have gradually brought to a critical pressure at shallower layers of the volcano to produce hydraulic fractures within the impermeable intervening caprock (Figures 6 and 9). Thus, increased injection rates and associated thermal effect may be linked to the opening of fractures to provide conductive medium for migrating hot water to reduce the mechanical strength of the fault, making it susceptible to failure (e.g., Figures 5 and 6). Similar mechanisms just related to pore pressure have been invoked also for earthquake triggering processes (e.g., Baccheschi et al., 2020;Chiodini, Cardellini, et al., 2020;Di Luccio, Ventura, et al., 2010;Miller et al., 2004).
We found that a minimum injection rate of 75-100 kg/s is required to cause seismic slip (Table 3 and Figure 5). The results suggest that the injection rate capable of causing the induced microseismicity at injection temperature of 350°C is near 100 kg/s (or ∼9,524 t/day). These observations confirm that the seismic unrest at CFc is explained by increased fluid injection rate from a deeper feeding system (e.g., Chiodini, Caliro, De Martino, et al., 2012;Chiodini, Paonita, et al., 2016;Chiodini, Todesco, et al., 2003;Coco et al., 2016;Rinaldi, Todesco, et al., 2010;Troiano et al., 2011). Similar to our work, Todesco, Rutqvist, Pruess, and Oldenburg (2003) and Todesco, Rutqvist, Chiodini, et al. (2004) increased the injection rates to simulate the HM and THM effects in relation to geochemical variation. According to Chiodini, Caliro, De Martino, et al. (2012), an injection fluid rate ∼30 times higher than the normal or background flow rate of the system (e.g., injection of 3,400 t/d or 39.7 kg/s of pure water or mixture of water and CO 2 ) is required to explain seismic unrests during 1983-1984, 1989, and 1995. Owing to the temperature capability of the TOUGH-REACT-FLAC 3D simulator, we could not adequately verify the contrast results of Afanasyev et al. (2015), where a fluid mass flow rate in the order of 50-100 kg/s at a depth of 5 km (and at ∼700°C temperature) and with an initial geothermal gradient of ∼120°C/km explains the ground deformation in the CFc area.
The results of this study indicated that seismicity can occur at a relatively greater depth mainly due to the hydromechanical effect with little or no thermal effect involved. This result is valid in regions with a high geothermal gradient such as CFc. The result may explain why microseismic events in the CFc region predominantly nucleated at depths lower than 3.0 km (Aster & Meyer, 1988;De Siena, Sammarco, et al., 2018;Di Luccio, Pino, et al., 2015). It is reasonable to interpret this result (depth-dependent THM effect) as due to the rheology of the caldera, specifically because of the deep brittle/ductile transition (∼3 km- Castaldo et al., 2019). This transition constraints the bottom of the seismogenic volume, allowing seismicity to migrate to shallower depths where thermal effects are enhanced (Chiodini, Paonita, et al., 2016;Chiodini, Vandemeulebrouck, et al., 2015).

Slip Analyses, Seismic Moments, and Magnitudes
The potential energy stored as elastic energy along fracture surfaces, such as faults, is dissipated during seismic slips and fault failures as other forms of energy, such as kinetic energy and heat energy. The moment magnitude (M s ) can be evaluated from the seismic moment (M 0 ) through the equation (Kanamori & Brodsky, 2001 is adequate for our work among the well-known empirical equations for computing the moment magnitudes, because it does not show saturation for large earthquakes (Hanks & Kanamori, 1979;Kanamori, 1977;Purcaru & Berckhemer, 1978).
The slip analyses reveal that displacements (slip distances) recorded along the fault during seismic events vary from millimeters to a few centimeters (e.g., Figures 5c, 7d, 8d, and 10e). Thus, low-to-moderate magnitudes (always M s < 3.0) are generated within the first 2.5 km depth. The distributions of the seismic slips along the fault (e.g., Figures 5c, 7d, and 10e) show that maximum slips are recorded near the middle of the fault, suggesting a link between the fault slip and evolution of the fault stress state and permeability from the fault's base as shown in Figure 11. This result is related to the rapid decrease in frictional strength on the fault's surface with slip followed by a gradual increase in frictional resistance which decelerates the slip motion, the observation which is similar to that in a slip-weakening process (Kanamori & Brodsky, 2001). This slip distribution has implications for spatial distributions of seismic events in the subsurface of CFc. In their discussion of the physics of earthquake ruptures and in analogy with this study, Kanamori and Brodsky (2001) noted that slip distributions can be obtained at various times while earthquakes are in progress, and slip can be heterogeneous both in space and time. Thus, analyzing slip distributions in this caldera can help gain further insights into its rupture dynamics especially when integrated with other monitoring data.
Increased injection rates and geological heterogeneities can affect the magnitudes of seismicity. Our results indicated that increasing injection rates could reduce required time to trigger fault slip, while enhancing the magnitude of seismicity (Figure 5c). Such episodic increased injection rates with potential increase in the magnitude of seismicity away from the background values have been attributed to intense magmatic degassing (Chiodini, Caliro, De Martino, et al., 2012;Chiodini, Todesco, et al., 2003), supply of hot-water-dominated fluids to shallow aquifer (Troiano et al., 2011), or feeding from a shallow magma source . Although the magnitudes of fault slip distance are higher under HM than THM in our results for Scenarios A and C (under varied injection rates and fault permeabilities, respectively), the thermal effect has more control on the magnitude of slip displacement in our Scenario D. The slip distances are more enhanced under THM at high caprock permeabilities (Supporting Information, Figure S5-compare panels (c) and (f)). The thermal effect also produced higher slip distance magnitude when the gap (temperature contrast) between the host rock temperature (assumed as 150°C) and injection fluid temperature (200°C -350°C) is significantly large (Figure 10e). Simulators traditionally require initial stability (equilibrium state) before gravitational loading (or hydromechanical and thermal loadings from volcano fluid injection as is the case in this work). Thus, the host rock formation temperature is based on the average surface temperature of 63.9°C at CFc (Vilardo et al., 2015) and an average geothermal gradient of ∼30°C/km, a typical value in a relatively stable nonvolcanic area within the first 3 km. It is remarkable that the highest injection temperature (350°C) with the largest temperature contrast gave the highest slip distance and vice versa, when the fluid injection rates are not more than 125 kg/s (Figure 10e). These results are consistent with the findings of Gan and Lei (2020) that high temperature contrasts led to the enhancement of fault slip distance. The enhanced slip displacement and corresponding magnitude of seismicity scaled with the highest shear stress drop (e.g., ∼0.53 MPa for highest injection temperature of 350°C-see Figure 10d). This relationship has been previously described for low-magnitude earthquakes (e.g., Drouet et al., 2011). However, we also found that slip distance began to diminish at very high injection rates (e.g., at 150 kg/s) when the injection temperature and thermal contrast are greater. Thus, our results seem to support the conclusion that the observed thermal effect at high injection rates might be due to fracture closure and fault healing at elevated temperatures (e.g., Cappa et al., 2019;McLaskey et al., 2012). We could not validate this claim because we solely considered time-dependent characteristics of fault friction and did not test chemical reactive effect (chemical component of the THMC model of Taron & Elsworth, 2009) in our numerical formulation. This is certainly the right direction for future works.
Our findings suggest that the dominant mechanism for seismic fault slip is hydromechanical in nature, but thermal effect may also play a role in affecting the timing of slip. While thermal effect could open up fractures faster and enhance permeability through thermally induced fracture aperture and shear dilation described in the preceding sections, it also causes earlier pressure dissipation which suppresses the pressure build-up required for slip distance enhancement ( Figure S5-compare panels (a) and (d), and panels (b) and (e)). The results also indicate that the caprock permeability plays a major role in determining the magnitude of seismicity. Besides prompting early fluid breakthrough and seismic fault failure, high caprock permeability also generated higher magnitude of seismicity under profound nonisothermal effect. The importance of caprock permeability cannot be overemphasized because it is the most significant factor influencing reservoir/caprock pore pressure distribution (Figures 6 and 9) and shear stress behavior at the injection reservoir rock interface (Zhou et al., 2015). The fault permeability and its contrast with the host rock matrix also determine the magnitude of seismicity. Our results suggested that low fault permeability and high-permeability contrasts (k m /k f in order of 10-1,000) produced higher magnitudes of seismicity (Figure 8d). Jasim et al. (2015) have also observed similar effects of k m /k f contrasts on the patterns of associated ground deformation at CFc. Thus, the major uncertainty in this work may come from our assumptions in the fault geometry as well as architecture and fault mechanical properties, which may be more complex in natural faults in volcanic regions (Cappa & Rutqvist, 2011). Future research must tackle a more accurate 3D description of the shallow lithological changes, for example, including tomographically described models that suggest the importance of lateral intrusions and vertical tectonic lineaments when modeling fluid propagation (De Siena, Sammarco, et al., 2018;Pepe et al., 2019;Siniscalchi et al., 2019). We have carefully summarized in Supporting Information ( Figure S5) the influence of each of the variables or factors discussed above on the fault displacements and respective magnitudes of seismicity based on our observations and deductions from the results of Scenarios A, C, and D to aid comparison between the HM and THM effects in this work.

The Onset of a Renewed Campi Flegrei Seismic Unrest-October 5, 2019
The fault slips (Figures 5c, 7d, 8d, and 10e) and depth-magnitudes (Figures 12a-12f) reported in this study, while comparable with the average magnitudes of the events recorded during the 1982-1984 magmatic unrest, can match neither their depth extension (down to 4.5 km) nor their maximum magnitudes (up to 3.8- Aster & Meyer, 1988;De Natale et al., 1991). In Figure 12, we compared spatial occurrences of the magnitudes of the events recorded by the Osservatorio Vesuviano from August 5, 2019, to December 5, 2019 (panel (a)) with those of the current modeled THM results (panels (b)-(f)). The comparison indicates that we can recover a realistic pattern of magnitude versus depths for the period August-November 2019, at the onset of the seismic unrest with the highest-magnitude earthquake swarms since the end of 1984 (December 6, 2019). Seismicity above M s = 1 nucleated between 0 and 1.5 km. The swarm on October 5, 2019, had its highest-magnitude earthquake (M s = 2.5) at a depth of 1.91 km, in agreement with most of the models that consider (1) high injection rates, (2) high caprock permeabilities, (3) low fault permeabilities, and (4) high thermal contrasts (maximum M s = 2.3 M w - Figures 12b-12e).
The results thus prove that the microseismicity above background seismicity (Chiodini, Selva, et al., 2017;Petrosino et al., 2018) at the caldera can be primarily modeled by hot-water injections in a capped reservoir, at least before December 6. This is only possible by including the unique properties of the Campi Flegrei caprock (Vanorio & Kanitpanyachaeron, 2015) into the model. The results show that depths and magnitudes can be better modeled in a hydromechanical (HM) approximation (Figure 12), in analogy with shortterm deformation signals (Coco et al., 2016). Even in the HM approximation, slip had to be modulated by (1) increased injection rates (panel (b)), (2) realistic geological heterogeneities and contrasts in fault and caprock permeabilities (panels (c) and (d)), and (3) relevant temperature contrasts between injected water and injection formation (panel (e)) to be realistic. until November 25, when the temperature started to recover. On the day of the December 6, 2019 swarm (M s = 3.1 M, depth = 2.33 km), the temperature had recovered to ∼112, within the standard trend of the last 2 years. The small deformation spike was stopped on the same date. While the depth of this earthquake is still comparable with those recovered in our model, such a high magnitude is well outside any simulation considering solely deep water injections. The temporal correlation with the rainfall activity producing the surficial thermal anomaly suggests a "double-injection model," from both a deep hot-water source and an injection of cold fluids from the surface. The deep subvertical structures of the faults under these craters (Bruno et al., 2017;Isaia, Vitale, Di Giuseppe, et al., 2015;Vitale & Isaia, 2014), recognized by geophysical imaging (Siniscalchi et al., 2019), are ideal to transform Pisciarelli into a natural injection well for cold rain waters.
Could such a double-injection model also reconstruct the latest high-magnitude earthquake swarm recorded at the caldera, which increased the maximum magnitude to 3.3M and maximum depth to 2.57 km-April 26, 2020, without the influence of shallow magma? If the situation was the same of 1983-1984, the December 6, 2019 swarm should have triggered the opening of the hydrothermal system, allowing an anomalous release of CO 2 associated with the earthquake (De Siena, Chiodini, et al., 2017). This CO 2 should come from the deep melt layer 7-8 km (Moretti et al., 2018;Zollo et al., 2008) or, especially, from a shallower magma batch near to the base of the model (3-4.5 km, Amoruso, Crescentini, Sabbetta, De Martino, et al., 2014;D'Auria, Pepe, et al., 2015). Including adequate increases of CO 2 concentration in the fluid injections is necessary to better answer this question. CO 2 emissions at Pisciarelli have dropped drastically just after the last swarm, on April 26, 2020 (Bollettino INGV-OV, May 5, 2020). Within 2 months, at the time of this paper submission, they had declined below the threshold set 11 years ago by the opening of the Pisciarelli crater. The behavior we described is typical of carbon capture storage under a highly impermeable caprock that rapidly recovers its integrity (Kampman et al., 2016), thanks to the unique ability of the caprock to withstand stress and recover from fracturing (Vanorio & Kanitpanyacharoen, 2015). Future modeling on the implications of permeability changes and chemical effects in the CFc caprock might thus be relevant to explain the evolution of the most important seismic unrest at the caldera in 35 years.

Conclusions
We show that TOUGHREACT can be coupled with FLAC 3D to explore the THM effects of hot-water injections in inducing instability of a volcanic caldera. This work has demonstrated that numerical simulations that take into account the complex geology of a volcano are not limited to modeling deformation and geochemical parameters: they can reconstruct seismic depths and magnitudes corresponding to an ongoing microseismic unrest. To do so, it was central to incorporate the material and mechanical properties of a caprock consistently inferred from rock physics and tomographic studies in our model parametrization. We hypothesized a two-step diffusion process in explaining the mechanism of caprock failure leading to fault instability. The first diffusion (fluid migration) between the underlying basal reservoir (injection formation) and the impervious caprock determines the failure of sealings, and the subsequent second diffusion between the fault and neighboring matrix determines the timing and magnitude of fault activation.
The initial low permeability of the caprock sealing allows pore fluid pressure to accumulate, which is critical to trigger seismicity by reducing effective normal stress state. The triggered caprock failure could lead to the breach of sealing, by enhancing caprock permeability at 5 orders of magnitude, through contrasting compressive and extensional forces at the basal and top parts of the caprock, respectively. Consequently, the uprising hot fluids could follow the created permeable channel to invade into fault.
High injection rates could prompt the onset of seismic slip with pronounced hydromechanical effects. Similarly, when rock cools down, the coupled THM influence enables the development of an early instability. However, the induced thermal effects could unload the compact state of the fault and open up fractures rapidly: this phenomenon leads to aperture and permeability enhancement accompanied by pressure dissipation, and, ultimately to a reduction in slip distance.
The magnitudes of seismicity generated are dependent on injection rates, thermal effect, structural heterogeneities from the caprock, and fault permeabilities. They are generally higher when fluid injection rates are increased and at lower fault permeabilities, but lower when caprock permeabilities are higher under isothermal conditions, due to less fluid invading into fault zones. The magnitudes and depths of the modeled microseismicity above the background threshold (1.0 M s ) are comparable to those recorded by seismic monitoring instruments at the start of the ongoing seismic unrest at Campi Flegrei. The comparison with the observed data suggests that processes analogous to natural carbon capture sequestration under an impermeable caprock could model the highest-magnitude seismicity (M s = 3.1 and M s = 3.3) of the unrest. This work has demonstrated the importance of numerical simulations in modeling seismo-volcanic unrests and quantifying the magnitudes of (micro)seismicity. The results of this study could be further constrained and integrated with relevant information from other techniques currently applied at the caldera, and then incorporated it into the routine volcano monitoring. These efforts shall undoubtedly improve short-to medium-term seismic and eruption forecasts, hazard evaluation, and mitigation actions beyond Campi Flegrei if integrated into a global database of volcanic unrests.

Data Availability Statement
The data sets for the results reported in this work have been archived in Zenodo and available at http://doi. org/10.5281/zenodo.4292032.