Modeling Intra‐Caldera Resurgence Settings: Laboratory Experiments With Application to the Los Humeros Volcanic Complex (Mexico)

Scaled analogue models explored the role of different boundary conditions in intra‐caldera resurgence processes. Models investigated the role of magma intrusion depth (ID) (Series 1), asymmetric and symmetric caldera collapse (Series 2), as well as the presence of existing discontinuities in the pre‐volcanic substratum (Series 3) on the style of caldera resurgence. Experimental results indicate that different IDs resulted in different resurgence styles and structural patterns, which evolved from piston‐like resurgence, for deeper intrusions, to intra‐caldera resurgent domes for shallower intrusions. Asymmetric collapse was typically accompanied by a tilted roof block above the emptied analogue magma reservoir, while inherited faults influenced significantly the deformation pattern of piston‐like resurgence. Experiments simulate many of the principal characteristics of calderas. We compare our modeling results primarily to the Los Potreros caldera nested within the Los Humeros Volcanic Complex, where the largest Mexico's Quaternary eruption occurred and which hosts an important geothermal field (eastern Trans‐Mexican Volcanic Belt). A structural field survey was conducted to identify the kinematics of faults within the caldera and outside the volcanic edifice. The Los Potreros caldera shows a sub‐orthogonal fault pattern strikingly similar to that of models deformed with shallow ID. We interpret this correlation as an evidence of similarity in dynamic processes, whereby modeling results would indicate a scaled ID of ∼4.5 km. The Acoculco caldera complex, in Mexico, shows a fault pattern similar to the Los Potreros caldera, and geological information corroborates the attribution of renewed magmatic pressure to similar IDs.

that has an obvious impact on caldera-forming systems (Cashman & Giordano, 2014;Kennedy et al., 2018). Styles of caldera collapse can be influenced by the degree of symmetry of collapse and deformation of the sinking roof block (e.g., trapdoor, piecemeal, downsag, and piston-like collapse), as well as the depth and dimensions of the magma chamber and the inherited structural setting (Acocella, 2007;Acocella et al., 2004;Cole et al., 2005;Lipman, 1997).
Post-caldera phases are often characterized by unrest and eruptions, as well as resurgence processes. Post-collapse resurgence is a process that induces uplift of a portion of the caldera floor, and is often associated with the input of new melts into shallow reservoirs (Galetto et al., 2017;Lipman, 1984;Marsh, 1984;Smith & Bailey, 1968). Intra-caldera resurgence characterizes many felsic calderas, but may also occur in mafic calderas (Cole et al., 2005;Lipman, 1984).
Caldera collapse and post-caldera resurgence are accommodated by the development of distinctive structural features (e.g., normal and reverse caldera ring faults, resurgent domes and associated apical normal faults), although some crucial aspects remain unclear or controversial. As for caldera collapse, unknowns and controversies regard the nucleation, propagation and dip direction of ring faults (see Geyer & Martí, 2014). Different models in fact have been proposed for internal fault configuration of caldera collapse, particularly regarding the dip direction of caldera-bounding (ring) faults (Geyer & Martí, 2014;Gudmundsson, 2007). The results of some experimental analogue models of collapsed calderas suggest that the sinking of brittle overburden overlying a magma reservoir is typically accommodated by early outward-dipping reverse faults and successive inward-dipping normal faults (e.g., Acocella, 2007;Geyer et al., 2006). In these models, the resulting ring fault system should therefore consist of concentric, peripheral normal faults delimiting more internal reverse faults, with the latter faults being often buried beneath the caldera products (e.g., Geyer & Martí, 2014;Geyer et al., 2006).
Post-caldera resurgence is accompanied by surface faulting, as well as seismic activity and the development of eruptive vents and sector collapse (e.g., Acocella et al., 2001;De Natale et al., 1997;Hill et al., 1985;Lipman et al., 2015;Orsi et al., 1999;Self et al., 1986;Tibaldi & Vezzoli, 2004). A recent model proposes that post-caldera resurgence is promoted by magma stagnation and renewed magmatic pressurization beneath a rheological barrier, which is made up of residual magma with viscosity higher than the underlying newly injected melts (Galetto et al., 2017). Intra-caldera uplift typically affects the central portion of collapsed volcanic systems, and is generally accommodated by normal faults developing at the apex of the uplifted zone (Bailey et al., 1976;Lipman, 1984;Self et al., 1986), with direction parallel to the long axis of the resurgent area (Galetto et al., 2017, and references therein), or resurgence faults may delimit the resurgent dome (e.g., Ischia and Pantelleria islands, Acocella et al., 2001).
Magmatic fluids stored in deep-seated reservoirs beneath caldera systems hold high potential for geothermal exploration and production (e.g., Fulignati et al., 1997;Giordano et al., 2014;Grigsby, 1982), especially in resurgent caldera settings that are often associated with the shallow emplacement of magma and the related development of hydrothermal systems within the caldera deposits (Lipman, 1984;Wohletz & Heiken, 1992). Circulation of geothermal fluids is typically controlled by the secondary permeability of faults connected to caldera and resurgence-forming processes (e.g., Garden et al., 2017;Giordano et al., 2014). Therefore, faults represent optimal targets for unlocking geothermal fluids (Martí et al., 2008). In addition, the characteristics of caldera resurgence is evidently depending on depth, size and position of the renewed magma source within the caldera system (Brothelande et al., 2016;Kennedy et al., 2012;Lindsay et al., 2001;Lipman, 1984).
Despite the knowledge achieved on this topic, some issues related to the deformation associated with caldera resurgence are still poorly known or controversial. These regard the architecture of both magmatic features and faults internal to the resurgent block and their relationships with emplacement rate and density of the resurgent magma, as well as with the existing caldera collapse faults. The latter can exert direct influence on the resurgence deformation outcome, but they in turn depend on a number of boundary conditions (i.e., asymmetry/symmetry of caldera collapse, presence of pre-existing faults in the pre-volcanic basement, etc.), which may complicate the study of the process. Besides an academic interest, the understanding of conditions controlling the structural architecture associated with caldera resurgence represents a research issue that bears primary interest for geothermal exploration and exploitation (Wohletz & Heiken, 1992).
In this study, we focus on the Los Humeros Volcanic Complex (LHVC) that is located in the eastern margin of the Trans-Mexican Volcanic Belt (TMVB), central Mexico. This system is characterized by two nested caldera collapses, the older is the larger Los Humeros caldera, which experienced the largest Quaternary eruption reported so far in the TMVB (Cavazos-Alvarez & Carrasco-Núñez, 2020), while the more recent is the inner and smaller Los Potreros caldera Ferriz & Mahood, 1984;Flores-Armenta & Gutiérrez-Negrín, 2011). The Los Potreros caldera hosts an important superhot geothermal system that is currently exploited for electric power production (Calcagno et al., 2018;Carrasco-Núñez et al., 2015). The Los Potreros caldera floor is affected by a peculiar pattern of faults made of interconnected fault segments collectively forming a sub-angular arrangement. In recent models, this setting has been attributed to intra-caldera resurgence and is inferred to control the circulation of geothermal fluids Norini et al., 2019Norini et al., , 2015.
Analogue modeling has been profitably used to investigate caldera resurgence processes (Acocella et al., 2000(Acocella et al., , 2001(Acocella et al., , 2004Brothelande & Merle, 2015;Galetto et al., 2017;Poppe et al., 2019;Urbani et al., 2020). On this basis, a new series of analogue models has been carried out to address the conditions that may lead to the structural pattern observed within the Los Potreros caldera. The present laboratory experiments offer an unprecedented procedure of modeling deformation that consisted of two distinct deformation stages given by an early caldera collapse stage and subsequent resurgence. The experiments have been conducted by systematically varying the depth of intrusion point, and testing different configurations of caldera resurgence, specifically piston-like and dike-sill processes, as well as the role of structures inherited from caldera collapse. The modeling results can shed light not only on the evolution and geometry of internal structures of the Los Potreros caldera, but can also provide some insight into other resurgent caldera systems worldwide. This information has obvious importance for reconstructing the caldera evolutionary model, as well as for identifying the internal fault architecture, which bears important implications for geothermal exploration and exploitation globally.
We begin to describe the results of the structural analysis carried on some faults affecting the Los Potreros caldera (LPC) floor. We then describe the results of an experimental series addressing the interrelation between the structural outcome and different initial boundary conditions (depth and shape of the new injected magma). Finally, we integrate the modeling results with structural observations to assess the tectono-magmatic controls on surface fault pattern of the LPC, as well as of other resurgent calderas in the region (particularly the Acoculco caldera complex in Mexico).

Geologic and Structural Setting of the Los Potreros Caldera
The LPC belongs to the LHVC, which is a prominent sub-alkaline, basaltic to trachytic and rhyolitic volcanic center (e.g., Carrasco-Núñez et al., 2017a;Lucci et al., 2020) located in the easternmost TMVB (Figures 1a  and 1b). The study area is characterized by a Palaeozoic crystalline basement grouped in the Teziutlan massif (Figure 1; e.g., Ferriz & Mahood, 1984). This crystalline basement is covered by a Mesozoic sedimentary succession deformed during the Late Cretaceous-Paleocene Laramide orogeny and now belonging to the Sierra Madre fold and thrust belt (e.g., Campos-Enríquez & Garduño-Monroy, 1987;Fitz-Diaz et al., 2018;Suter, 1987). Oligocene to Miocene magmatism is mainly characterized by scattered distributed felsic intrusive bodies and by the dominantly andesitic magmatism of the Cerro Grande volcanic complex (Gómez-Tuena & Carrasco-Núñez, 2000) and the Cuyoaco and Alseseca andesites (see Lucci et al., 2020, and references therein). Late Pliocene to Early Pleistocene magmatic activity is represented in the study area by the voluminous Teziutlan andesites. These Pliocene andesites correlate with the thick andesitic layers of the subsurface stratigraphy of the Pleistocene-Holocene LHVC (Carrasco-Núñez et al., 2017b), which experienced at least two large caldera collapse events (Ferriz & Mahood, 1984). The older caldera collapse produced the 18 × 15 km wide Los Humeros caldera associated with the emplacement of the Xaltipan ignimbrite, whose dense rock equivalent (DRE) volume was first estimated in 115 km 3 (Ferriz & Mahood, 1984) and more recently revised in ca. 285 km 3 (Cavazos & Carrasco-Núñez, 2020) ( Figure 1b). Therefore, it represents the largest Quaternary eruption in the TMVB. The second caldera collapse event led to the formation of the bimodal andesitic-rhyodacitic Zaragoza ignimbrite (15 km 3 DRE) as well as the nested 9 × 10 km-wide LPC, which developed close to the south-western margin of the Los Humeros caldera   Norini et al. (2019), andC. Geológico-Minera (2002)]. Sites of fault slip-data measurement are indicated. The lower hemisphere stereoplots report the stereographic projection of faults (thin black lines), magmatic dikes (red lines), joints (thin blue lines), axial planes of folds (thin green lines), and fold axes (green circles); the small black arrows indicate the movement direction of the fault hanging wall for each slickenline measurement, while the small white arrows represent kinematic indicators in which a reliable sense of shear was not observed. Divergent and convergent open arrows indicate respectively the σ 3 and σ 1 direction computed through the right dihedron inversion method (Angelier, 1979) (structural data in Table 1). Stereonets of Sites 1-7 are reported in Figure 2. Site 14 is located east of the mapped area. (c) Geological section across the LHVC (redrawn after Norini et al., 2015). Dashed black lines, Los Humeros caldera bounding fault, blue line, Los Potreros caldera bounding fault; red line, intracaldera Los Potreros fault. structures thus affect the early Teziutlán volcanic succession, which is also overlain by the Xaltipan and Zaragoza ignimbrite deposits (Carrasco-Núñez et al., 2017a;Figure 1c). Radiometric dating of caldera collapse events has been recently constrained using two geochronometers (zircon U/Th and plagioclase 40 Ar/ 39 Ar dating), yielding ages of 164 ka and 69 ka for Xaltipan and Zaragoza ignimbrites, respectively . The Pliocene(?)-Pleistocene Teziutlán andesites are instead dated ( 40 Ar/ 39 Ar) at 1.46-2.61 Ma (Carrasco-Núñez et al., 2017a). Interpretative geological cross sections show both caldera structures as being controlled by either inward-dipping (Cedillo-Rodríguez, 1999;Ferriz & Mahood, 1984), or outward-dipping faults (Calcagno et al., 2018;Carrasco-Núñez et al., 2017b;Norini et al., 2019Norini et al., , 2015. In map-view, the south-eastern margin of the Los Humeros caldera is bounded by a ca. NE-trending linear scarp referred to as "Los Humeros scarp," whereas its south-western rim is conceivably controlled by a NW-trending nearly rectilinear fault that is hinted from the clear alignment of several volcanic vents (Carrasco-Núñez et al., 2017a;Figures 1b and 2). This setting may suggest that the geometry of the Los Humeros caldera rim was controlled by reactivation of inherited regional structures during the collapse (Norini et al., 2019;Liotta & WP4 Working Group, 2019).
The LPC is delimited by the "Los Potreros scarp" and the "Oyameles scarp" along its eastern and western margins, respectively (Norini et al., 2015; Figure 2). Such exposed caldera collapse scarps vanish abruptly laterally, or are buried beneath the thick syn-to post-caldera products. The LPC's floor is affected by a system BONINI ET AL.   Urbani et al. (2020). Caldera rims are from Norini et al. (2015). Sites of fault slip-data measurement are indicated. Symbols are as in Figure 1. Coordinates of the sites of structural measurements are shown in Table 1 of active normal faults, and subordinately reverse faults, defining a main NNW-trend from which curved fault segments depart with a trend varying from NE to E-W (Norini et al., 2015). Importantly, these faults represent the main target of geothermal exploration and production, and are inferred to delimit intra-caldera sectors with similar vertical uplift (Norini et al., 2015). Differential uplift among these blocks has been related to intra-caldera magma resurgence (Norini et al., 2019(Norini et al., , 2015Figure 2), which is also suggested by the presence of subsurface rhyolitic domes that have been drilled by geothermal wells (see geological section in Carrasco-Núñez et al. (2017a).
Some studies suggest that the magma chamber-that is expected to source magma resurgence-is currently partly solidified, with present-day temperatures of the order of 600°C-650°C (5 to 7-8 km depth; e.g., Verma, 1983Verma, , 2000. However, a recent study suggests that the post-caldera magmatic plumbing system would instead consist of a complex arrangement of multiple magma storage layers located at different levels between a deep basaltic reservoir (ca. 30 km) and shallow-crust conditions (depth of ca. 3 km) . In this scenario, recent deformation and uplift within the LPC is related to the ascent of several, very shallow (<1 km) multiple magma bodies (producing lava domes and cryptodomes), which deform small-scale areas at the surface (areal extent ∼1 km 2 ) showing normal faulting at their top .

Methods
Tectonic elements (fault and fracture sets, fold axes) and magmatic features (dikes) have been measured within the LHVC, as well as in the pre-volcanic substratum outside the LHVC (Table 1). Structural field analysis focused on the LPC's floor, which is the most active part of the volcanic edifice that seemingly records the recent/ongoing tectono-magmatic deformation. Kinematic indicators in the LPC were measured from major intra-caldera fault segments, as well as on associated minor faults. Fault slip data (n = 140; data available at http://gemex.igg.cnr.it/documents/231) contain fault orientation (fault dip and dip direction), and sense of fault motion (rake) obtained from fault kinematic indicators (Angelier, 1994;Doblas, 1998). Fault-slip data have been subdivided into those measured on (i) the major fault plane (and fault core), (ii) in the fault damage zone, and (iii) on other mesoscopic faults distant from a major fault (Table 1). Kinematic analysis of a fault population was performed through the graphical right dihedron angle method of Angelier and Mechler (1977) and Angelier (1979), which allow the identification of the areas in a stereoplot where the principal stress axes (σ 1 , σ 2 , and σ 3 ) have the maximum probability to lie. This method was integrated with the condition of orthogonality among the stress axes (Caputo & Caputo, 1988).

LPC
Evidence for recent faulting affecting the volcanic edifice has been dominantly found within the LPC. In particular, the central part of the LPC is affected by a system of intra-caldera faults that are characterized by fresh fault scarps defining a sub-radial pattern. More specifically, such a fault system exhibits a complex pattern with fault segments showing two main orientations, particularly (1) faults trending NNW to N-S and dipping west to southwest (Maxtaloya and Los Humeros faults, as well as some other minor faults) and (2) faults striking from E-W to NE-SW (Las Papas and Arroyo Grande faults; Figure 2). In particular, the minor E-W to NE-SW-trending faults branch-off at high angle from the main NW-striking Maxtaloya and Los Humeros faults, showing a curvilinear shape and defining a sort of "horse-tail pattern" (Figure 2

Note.
Letters "a-c" beside a site number indicate the inferred chronological cross cutting relationships, with "a" being the older event. See text for details.  (Angelier, 1979) direction (azimuth 246°) that trends sub-orthogonal to the local strike of this structure (Figures 2, 3a, and 3b and  Table 1). Dislocation of reference markers allowed us to determine the dominant type of movement (extensional or reverse), but reliable fault kinematic indicators were not observed, thereby preventing application of stress inversion to these fault populations. Furthermore, two sites of outcrop-scale faults measured in the damage zone of the Arroyo Grande fault show a dominant N-S trend (Sites 7a and 7b, Figure 2). Kinematic indicators collected on Site 7a reveal dominant pure dip-slip normal faulting and an ENE-oriented extension direction (azimuth 252°) that trends sub-orthogonal to the local strike of the Arroyo Grande fault ( Figure 2 and Table 1). Both normal and reverse faults were observed at Sites 4, 6, and 7, although their chronological relationships are unclear ( Figure 2 and Table 1). Reverse fault occurrence would reflect the curvature of steep resurgent-related fault planes (e.g., the Arroyo Grande fault; Norini et al., 2015Norini et al., , 2019.

Pre-volcanic Substratum
The pre-volcanic substratum is mostly made of Mesozoic limestone and shale that have been folded and faulted during the Laramide orogenesis. Structural measurements carried out at selected sites in the pre-volcanic substratum evidence that Laramide-related deformation consists of ca. NW-SE-trending thrust folds and associated joint systems, of which cross joints sub-orthogonal to fold axis are often dominant (Site 8; e.g., Norini et al., 2019). Application of stress inversion produces consistent maximum horizontal stress axes trending between NNE and ENE (Sites 8, 10, 11a; Figure 1b and Table 1).
The pre-volcanic substratum also records distinct deformational events superposed onto the early thrust and fold structures. In particular, the pre-volcanic basement is affected by NNE to ENE-trending faults showing different fault kinematics, namely strike-slip, transpressive, and normal-oblique-slip (Sites 9,13,15,16,17,18; Figures 1b and Table 1). At some places, transverse strike-slip to transpressive faults are superimposed onto thrust faults, and the normal-oblique-slip faults post-date the strike-slip to transpressive faults. This complete chronology of deformation has been observed at the thrust anticline exposed west of Tenextepec (Site 11; Figure 1b). There, the three fault populations show paleo-stress axes with similar orientation, possibly suggesting the occurrence of multiple tectonic permutations (Table 1). At the same locality, ca. N-S-trending normal faults displace cemented slope breccias (Site 12), probably revealing a late normal faulting event ( Figure 2 and Table 1). The pre-volcanic basement was also deformed by ca. NW-SE-trending faults (Sites 14-16; Figure 1b). In some cases, the NNE/ENE-and NW-SW-trending faults have been found to displace each other (Site 16), suggesting that they developed during a same deformation stage. This surface deformation may reflect the reactivation of inherited, subsurface NW-SE and NE-SW-trending fault sets that have been identified in this region (García-Palomo et al., 2018;Liotta & WP4 Working Group, 2019). This hypothesis is further supported by the observation that rocks of pre-volcanic substratum are intruded by NE-SW and NW-SE-trending magmatic dikes (Sites 9, 10, 15, 17; Figures 1b and Table 1). Magmatic dikes have been found to exploit NE-and NW-trending discontinuities, and some volcanic centers are frequently aligned along NE directions, suggesting a control of deep-seated (Laramide) structures with similar orientation in the pre-volcanic substratum (Norini et al., 2019).

Modeling Strategy
Analogue modeling was carried out at the Tectonic Modeling Laboratory of the CNR-IGG and of the Department of Earth Sciences in Firenze, Italy (raw data of models presented in this study are freely available to download from Maestrelli et al. (2020) at the link http://doi.org/10.5281/zenodo.3882130). Models were appropriately scaled to nature (see Appendix A for the complete scaling procedure) and designed to test the influence of both depth and modality of renewed magma intrusion on the deformation observed in map view and in cross section of resurgent caldera systems. Our experimental series has considered a two-phase deformation history consisting of (1) initial caldera collapse (reference model SC-01,  The second phase of intra-caldera magma resurgence has considered the cases when renewed magma injection occurs into the early magma reservoir that generated the collapse (De Silva et al., 2015;Kennedy et al., 2012;Marsh, 1984;Smith & Bailey, 1968), or into shallower intrusions (Fridrich et al., 1991;Hildreth, 2017;Kennedy et al., 2016) (Figure 4). For this reason, in some experiments the depth of intrusion point was varied systematically and positioned at progressively shallower depths in order to investigate the role of magmatic intrusions on intra-caldera resurgence (Series 1). In other models, the analogue magma was injected into the same reservoir that was previously evacuated to produce caldera collapse, according to a piston-like collapse and resurgence mechanism (Series 2). Models have been generally deformed applying a symmetric collapse, but some models have also explored the role of asymmetric (trapdoor) collapse that produces structures that may be potentially reactivated during caldera resurgence (Series 2). Furthermore, some models have considered the role of sub-orthogonal, inherited discontinuities in the pre-volcanic substratum (see Section 3.2.2) in the development of caldera collapse and resurgence (Series 3).

Model Construction, Deformation, and Materials
The setups of models presented in this study are summarized in Figures 4c-4e and Table 2. The model setup consisted of a brittle overburden overlying a low-viscosity fluid simulating magma in nature. The brittle overburden had a standard thickness of 6 cm, and in few models 4 or 8 cm. The low-viscosity fluid was confined within a 2 cm-thick cylindrical container made of plaster with a diameter of 12 cm, which mimics a magma chamber in nature ( Figure 4a; Series 1 and Series 2). In few models, the analogue magma chamber was delimited by two straight margins simulating pre-existing orthogonal discontinuities in the pre-volcanic substratum (Section 4.1; Figure 4b; Series 3). Such discontinuities were artificially introduced also into the brittle overburden as vertical cuts made with a knife to reproduce the presence of pre-existing structures in the pre-volcanic substratum, whereby they stop at 2 cm depth (e.g., Maestrelli et al., 2021; Figure 4b).
Models were built in a Plexiglas box with internal dimensions 33 × 38 × 15 cm ( Figure 4).

Table 2 Parameters of Analogue Models Discussed in This Study
The brittle overburden was simulated by a mixture of Quartz and Feldspar sand (proportions 70:30 in weight), a Mohr-Coulomb granular material with density of 1,408 kg m −3 , angle of internal friction of ∼40° (coefficient of friction μ ∼ 0.83) and cohesion of ∼10 Pa (Del Ventisette et al., 2019;Montanari et al., 2017a). Sand mixture was poured in the model and then shaved with a leveling apparatus to produce layers of colored sand used as passive markers to visualize internal model deformation. The low-viscosity fluid is the pure vegetable polyglycerine-3 (PG3) manufactured by Spiga Nord S.p.A. This material has a Newtonian behavior, a density of 1,190 kg m −3 and viscosity of 17 Pa s (Montanari, et al., 2017b).
The two-phase deformation started with the evacuation of the analogue magma through a cylindrical PVC pipe (9 mm in internal diameter) located below the center of the analogue magma reservoir, a process that induced the collapse of a caldera (phase 1). The fluid was let flow out under the weight of the brittle overburden, so that the outflow velocity of the experimental magma was similar in all models. The process was stopped when only a thin layer of viscous material remained in the model reservoir (approximately the 65% of initial analogue magma). Afterward, intra-caldera resurgence (phase 2) was simulated by re-intruding the analogue magma (from the center of the caldera) at different depths into the brittle overburden (Series 1 and partly Series 3), or into the fluid container evacuated during caldera collapse (Series 2 and partly Series 3) (Figure 4). In particular, four models tested the effect of different injection depth (Series 1), which was systematically varied from 3 to 5 cm depth over an overburden thickness of 6 cm (Figures 4c-4e). In these models, the fluid was re-intruded into the brittle overburden through a 4 mm diameter pipe, which was BONINI ET AL.
10.1029/2020JB020438 11 of 32  Table 2). The latter models incorporate artificial discontinuities that extend up to 2 cm depth into the brittle overburden. Sketch showing map-view and lateral view of experimental setups of models of (c) Series 1, (d) Series 2, and (e) Series 3. The red dot indicates the inlet position both in map view and cross-section. ID, injection depth.
inserted into the largest tube used to evacuate the material during the previous caldera collapse stage. Injection of fluid simulating magma was driven by a piston controlled by an electric stepping motor predominantly at a rate of 80 cm h −1 (Table 2). Five models tested the deformation patterns produced by a piston-like process (Series 2), in which the fluid simulating magma was intruded into the previously drained reservoir through the same tube used to remove the experimental magma (Figures 4c-4e and Table 2). In the models that explored the role of asymmetric collapse (Table 2), the latter was induced by tilting slightly the top of the analogue magma chamber, so that the non-uniform overburden thickness generated a differential vertical loading over the surface of the analogue magma (lower schematic cross-section in Figure 4d). Intrusion of experimental magma was carried out at rates of 80 and 40 cm h −1 depending on the model (Table 2). Two models explored the role of two pre-existing, orthogonal discontinuities in the pre-volcanic substratum that delimited the analogue magma chamber and extended into the overlying brittle sand pack (Series 3). Models also incorporated a vertical discontinuity in the brittle overburden centered on the injection point (Figures 4c-4e and Table 2).
Some models were dedicated to test model setups and reproducibility of results. After deformation, models were covered by dry sand to protect the final surface topography. Models were then impregnated in water and frozen to allow sectioning without disturbing the model. We report below only those that are used to discuss the final results of the modeling.

Monitoring and Analysis of Model Deformation
Deformation of models was monitored by automatically taking high-resolution top-view photographs (12 megapixel) at constant time intervals of 120 s. Digital elevation models (DEMs) of models surface were elaborated at the end of both caldera collapse and resurgence phases through Agisoft Photoscan® software, which is based on photogrammetric techniques that are commonly used for monitoring geological laboratory experiments (e.g., Donnadieu et al., 2003). In particular, dedicated high-resolution photographs (18 megapixel) were acquired around the target model from fixed viewpoints. This procedure allowed to obtain a three-dimensional rendering (point dense cloud) of model surface and the interpolated DEM. DEM resolution ranges from a minimum of 0.356 mm/pix to a maximum of 0.152 mm/pix, with an average resolution of 0.251 mm/pix. Markers were placed at fixed positions on the model setup, allowing the equal scaling of all DEMs. The zero-reference level for elevation was arbitrarily attributed to flat, undeformed model surface. Beside a sharper view of fault scarps, the geometry attributable to resurgence can be efficaciously visualized by subtracting the DEM at the end of resurgence phase from that at the end of caldera collapse stage.
Deformation of model surface was also implemented with Digital Particle Image Velocimetry (DPIV) analysis through the PIVlab algorithms developed by Thielicke and Stamhuis (2014) for MATLAB®. DPIV is a non-invasive, quantitative and qualitative technique commonly used for visualization of particle flow. In particular, PIVlab allows the calculation of particle displacement vectors over model surface between two selected photographic frames, and allows the quantification of two-dimensional model deformation.

Simplifications and Limitations of the Modeling
Analogue models necessarily simplify the complexity of natural geological systems owing to the unfeasibility to reproduce all the components involved in the process. This difficulty is particularly evident in magmatic systems, where temperature-dependent processes are important but unfortunately cannot be addressed. This implies the application of important limitations to the modeling, which basically regard (i) the rheological properties of the analogue magma, (ii) the role of volatiles and (iii) the variations in rheology during cooling and crystallization of magmas, and (iv) the rates of magma emplacement (e.g., Corti et al., 2003;Merle, 2015). Our modeling procedure ignored thermal effects as magma was simulated by a low-viscous fluid, whereby the complex thermo-mechanical behavior of the system was approximated to a purely mechanical deformation. In addition, the boundaries of the experimental magma chamber were rigid (undeformable).
Another shortcoming regards the experimental fluid injection apparatus, which consists of a 4 mm pipe scaling down to 0.4 km in nature. Natural magma conduits have a diameter generally less than 50 m, and exceptionally can exceed 100 m (e.g., Aravena et al., 2018), so that our injection apparatus can only simulate a zone of pervasive magma transport (e.g., Benn et al., 2000).

Description and Discussion of Model Results
The early caldera collapse in the models was initially accommodated by the development of reverse faults dipping outward away from the ensuing depression. The development of the reverse fault system was not synchronous but, after initial down-sagging, deformation nucleated at a specific point and from there propagated laterally following a circular trajectory. With increasing collapse, deformation became accommodated by a system of normal faults dipping toward the caldera center (see below). These faults widened the depression and created an external ring fault system ( Figures 5-7, end of collapse). This evolution is basically similar to what is described in previous analogue models of caldera collapse (e.g., Acocella, 2007;Geyer et al., 2006;Holohan et al., 2013;Martí, et al, 1994;Roche et al., 2000). The presence of straight discontinuities bounding the analogue magma chamber influenced only partly this evolution, which will be explained in Section 5.3. The experimental results are described below starting from Series 1, which addresses the influence of intrusion depth (ID) on the deformation outcome of model resurgence. We then describe the experiments exploring piston-like resurgence ensuing symmetric and asymmetric caldera collapse (Series 2), and finally the models built with pre-existing discontinuities (Series 3).

Variation of Intrusion Depth (Series 1)
In Model LR-01, the fluid simulating magma was intruded at 5 cm depth and lifted the whole caldera floor (see model top-views at the end of caldera collapse and resurgence; Figure 5a). Lifting of caldera floor is illustrated in Figure 5a (central panel) as the topographic difference between the DEMs at the end of resurgence and at the end of caldera collapse. Analysis of model cross section reveals that deformation was accommodated by a deep sill intrusion that essentially behaved as a piston-like system, which reactivated the peripheral caldera normal faults as reverse faults (Figure 5a).
In Model LR-02, the experimental magma was injected at 4.5 cm depth (Figure 5b). At surface, the experimental caldera floor was deformed by a resurgent dome affected by a system of normal faults trending at high-angle from each other. In comparison to the former model, the resurgent area was smaller and more asymmetric. Interestingly, the model DEM shows that resurgence-related faults delimit sectors with different elevation (central panel in Figure 5b). In particular, the resurgent dome is characterized by a depressed, fault-bounded apical sector, which is in turn delimited by sub-orthogonal fault segments dipping in the opposite direction to the lowered sector. The model cross section reveals the occurrence of a thick saucer-shaped sill that has moderately to steeply dipping inclined sheets at its margins, from which inward-dipping reverse faults nucleate (Figure 5b). Resurgence of analogue magma is thus accommodated by a block forced upwards by steep, inward-dipping reverse faults. The section also shows that the outward-dipping normal fault, orthogonal to the apical depression, dies out at depth above one of these conjugated reverse faults (similar to deformation features associated with forced folding in Montanari et al. (2017b) and Poppe et al. (2019).
Model LR-03 was deformed with a shallower intrusion point, settled at 4 cm depth beneath the caldera floor ( Figure 5c). In top-view, the experimental resurgent dome was deformed by a few normal fault segments, and the analogue magma extruded at surface near the dome rim. The resurgent dome was smaller than that developed in Model LR-02. The topographic difference illustrates a clear bulging of the caldera floor, with a secondary inflation on the left-hand side that is related to the extrusion of the analogue magma (central panel in Figure 5c). The model cross section shows that fluid intrusion resulted in a small sill above the fluid inlet, from which departed two conjugated, steeply dipping sheets. One of them was more developed, and produced the surface extrusion of the experimental magma. The conjugated sheets are associated with reverse faults delimiting a lifted block that is narrower than that developed in Model LR-02.
Finally, in Model LR-04 the analogue magma was intruded at 3 cm depth (Figure 5d). At surface, the caldera floor was deformed by a small resurgent dome affected by a diffused deformation centered over the dome and consisting of small, closely spaced radial fractures and faults. Like Model LR-03, the experimental mag-BONINI ET AL.

10.1029/2020JB020438
14 of 32 ma pierced the external carapace of the dome and reached the surface. The differential topography shows a localized bulging of the caldera floor (central panel in Figure 5d), which is the smallest surface among those of this series of models. The model cross-section discloses that fluid injection produced a sub-vertical dikelike intrusion. This intrusion was delimited by two conjugated reverse faults that also created the surface resurgent dome.

Role of Symmetric and Asymmetric Caldera Collapse on Piston-Like Resurgence (Series 2)
Model PR-01 consisted of an initial phase symmetric caldera collapse, which was accomplished by the typical development of early reverse faults and subsequent peripheral normal faults (Figure 6a). Injection of the analogue magma into the reservoir lifted the whole caldera floor, as also highlighted by the topographic difference between final and initial stages of caldera resurgence (central panel in Figure 6a). Interpretation of the top-view photograph also suggests that the external normal faults have been reactivated as reverse BONINI ET AL.

10.1029/2020JB020438
15 of 32 Models PR-02 to PR-05 investigated the role of asymmetric caldera collapse. Model PR-02 experienced fluid injection at a rate of 80 cm h −1 (Figure 6b), while Models PR-03 to PR-05 were deformed with an injection velocity of 40 cm h −1 and explored the role of different overburden thickness (Figure 7 and Table 2). The collapse phase resulted in the typical development of reverse and normal faults, but this process was accompanied by a marked tilting of the roof block, with geometry similar to trapdoor collapse (Figures 6b and  7). However, asymmetry was less pronounced in model PR-05, presumably because differences in vertical stress at the base of the experimental roof block were less important than in models built with a thinner overburden (Models PR-03 and PR-04). Tilting of the experimental trapdoor block was accommodated by an array of normal faults (mostly in Models PR-03 and PR-04) that bounded its lifted margin and dipped oppositely to the direction of block tilting. These faults developed near to the caldera rim, dipping in the opposite direction to the ring-normal faults.
Model top views at the end of the resurgence phase show that the peripheral normal faults inverted their movement and were reactivated as reverse faults, similarly to Model PR-01. On the contrary, the normal fault that hinged the trapdoor block did not reverse its movement but increased its vertical throw. Small newly formed fault segments developed in the center of the trapdoor block of some models (e.g., Model PR-03; Figure 7). The experimental caldera floor was thus affected by a number of fault-bounded blocks showing differential uplift, which is highest on the lifted margin of the trapdoor block. This setting is well BONINI ET AL.

10.1029/2020JB020438
16 of 32 outlined by the differential topography between caldera collapse and resurgence, where the normal fault visibly delimits a more uplifted sector that corresponds to its footwall (central panel in Figure 6b).
The model cross sections reveal the asymmetric inflation of the injected analogue magma below the tilted trapdoor block, and a similar internal fault architecture (Figures 6b and 7). Also models addressing the variation in overburden thickness (models PR-03 to PR-05) produced basically similar deformation patterns (Figure 7). The cross sections support the interpretation that the external normal faults reversed their movement to accommodate the lifting of the caldera floor and trapdoor block. However, most of deformation was taken up by the faults bounding the apical part of this block, which gained further vertical displacement.

Role of Pre-existing Discontinuities on Resurgence Deformation Patterns (Series 3)
The occurrence of linear discontinuities in the pre-volcanic substratum (see Section 4.2) affected caldera collapse at some extent. Whereas the typical progression from early reverse faults to late normal faults developed in the circular part of the chamber, the rectilinear segments prevented the reverse fault development. The reason is that the pre-existing vertical discontinuities were reactivated since the very early stages of collapse, accommodating most of the deformation. We present the results of two models built with this configuration that have tested different injection depths of the experimental magma.
Model DR-01 tested the case in which injection was applied to the base of the brittle sand overburden (6 cm depth). As described earlier, the initial phase of symmetric caldera collapse was completed through the development of rectilinear caldera rims, and the usual progression from reverse to normal faulting in the circular sectors (Figure 6c). Injection of the analogue magma lifted the whole caldera floor, a process that was achieved through the reactivation of all the peripheral normal faults. A symmetric uplift of the caldera floor is also demonstrated by the differential topography of the caldera floor (central panel in Figure 6c). The model cross section validates the hypothesis of piston-like resurgence of the caldera floor, as well as the consistent reactivation of the outer caldera normal faults as reverse faults (Figure 6c). Instead, we observed no reactivation of central discontinuity.
Model DR-02 investigated the effect of a shallower intrusion point, which was settled at 4.5 cm depth beneath the caldera floor ( Figure 6d). In contrast with the previous model, the experimental resurgent dome affected a relatively small portion of the caldera floor. The apical part of the dome was affected by a diffused deformation, and by the extrusion of the magma-simulating fluid near the forced fold rim. Although only few short fault segments could be detected, these show a clear radial distribution. The differential topography of caldera floor shows a clear bulging of the caldera floor and a secondary inflation on the right-hand side that is associated with the surface extrusion of the magma analogue (central panel in Figure 6d).
The model cross section reveals that fluid intrusion created two steep, conjugated dikes that developed above the fluid inlet. One of them reached the surface and extruded the fluid, but this occurrence is not captured by this section. The cross-section structural pattern is similar to that of Model LR-03 (Series 1), where the dikes are associated with reverse faults delimiting a portion of the experimental caldera that is forced upwards (cf. Figures 5c and 6d). However, Model DR-02 experienced a slightly deeper ID (4.5 cm) than Model LR-03 (4 cm). Injection depth of Model DR-02 was equal to that of Model LR-02, but the final model outcomes were rather different, Model LR-02 being characterized by a saucer-shaped sill (cf. Figures 5c and 6d).

Results of DPIV Analysis
DPIV analysis was applied to all model duration, but it is reported only for time intervals conveniently showing the deformation outcome of each model (Figure 8). It is worth recalling that DPIV outlines the horizontal component of displacement over the model surface. DPIV results outlines effectively the different evolutionary models observed for Series 1. In particular, the velocity field shows a rather homogenous radial pattern that is consistent with the piston-like resurgence characterizing Model LR-01 ( Figure 8).
DPIV results match well also the deformation patterns resulting from the shallowing of ID. In particular, the distribution of velocity vectors highlights the occurrence of different fault-bounded blocks in Model LR-02, as well as the localization of deformation in models LR-03 and LR-04 ( Figure 8). The latter models show again a radial, centrifugal velocity pattern, in which there is an area of increased velocity that corresponds to the region of future extrusion of the fluid simulating magma.
Regarding Series 2, DPIV results show the occurrence of a radial pattern for Model PR-01, which is not, however, centered on the intrusion, possibly because of the slight tilting of the roof block. Model PR-02 shows instead a more asymmetric velocity pattern (Figure 8). These results are basically consistent with the symmetric and asymmetric collapse and resurgence imposed to Model PR-01 and Model PR-02, respectively (see Figures 6a and 6b).
Series 3 illustrates the control on deformation pattern exerted by both pre-existing discontinuities and ID. DPIV analysis of Model DR-01 shows the typical radial velocity field produced by resurgence, but velocity distribution is asymmetric in that horizontal velocities are highest in correspondence of both the pre-existing faults and on the opposite side ( Figure 8). The velocity pattern obtained for Model DR-02 is instead similar to those resulting from shallow ID (i.e., models LR-03 and LR-04; Figure 8).
A feature common to many models is that horizontal velocity shows lower magnitudes in the apical part of the experimental resurgent dome, or in the central part of the resurgent piston-like caldera floor (where velocity has a dominant vertical component that is not captured by DPIV analysis). Horizontal velocity is instead maximum at the external regions of resurgent domes, most likely in association with the outwards movement of the inward-dipping reverse faults that delimit the dome or lift the caldera floor.
BONINI ET AL.

Role of Injection Depth in the Tectono-Magmatic Deformation
Our laboratory models simulated the emplacement from different depths of magmatic intrusions into collapsed calderas. The resulting experimental intrusion geometries are similar to those obtained in previous analogue models exploring magma emplacement at shallow crustal levels (Galland et al., 2007(Galland et al., , 2014(Galland et al., , 2009Mathieu et al., 2008;Montanari et al., 2017b;Montanari et al., 2010b;Montanari et al., 2020;Poppe et al., 2019;Román-Berdiel et al., 1995). Previous analogue models of caldera resurgence systems investigated only the resurgent phase (Acocella et al., 2000(Acocella et al., , 2001Urbani et al., 2020), or inflation-deflation cycles produced by a balloon embedded into brittle materials (Troll et al., 2002;Walter & Troll, 2001). Our models are instead the first to investigate the development of resurgent-related tectono-magmatic structures within a caldera structure that formed under different boundary conditions (e.g., circular magma chamber, presence of rectilinear discontinuities, symmetric and asymmetric collapse).
The shape of experimental magmatic intrusions was strongly influenced by the different injection depths, therefore by the differential loading exerted by the brittle overburden on the analogue magma. In order to illustrate this behavior, the ID has been contrasted with the resurgence-induced deformation area (RIDA) observed on the model caldera floor (Figure 9a). Models show that the resurgent area on the experimental caldera floor clearly decreases with decreasing injection depth. This variation is illustrated in the graph ID versus RIDA of Figure 9b. The variation of resurgent area on the caldera floor is matched by the development of different intrusion geometries at depth, which are sketched in Figure 9c. For deep injection, the analogue magma formed laccoliths or saucer-shaped sills, which inflated and forced upwards the whole or the great part of the caldera floor. For shallower intrusion, deformation of the overlying brittle overburden consisted of steep dike systems whose superficial effect were more limited as only part of the caldera floor was affected by doming. Another important observation is that models incorporating a deep intrusion of fluid produced piston-like geometries that lifted uniformly the caldera block (Figures 5a, and 6a-6c). The steep dike systems associated with shallow resurgence of experimental magma was instead accompanied by the development of two conjugated reverse faults that folded and lifted the interposed portion of the caldera block (Figures 5b-5d and 6d). Noteworthy, folding of the roof block becomes more pronounced with decreasing ID.
In natural systems, "forced folding" is a primary mechanism creating space for magma emplacement at shallow crustal levels, and it is normally associated with the development of brittle deformation in the overlying overburden (e.g., Hansen & Cartwright, 2006;Jackson et al., 2013;Pollard & Johnson, 1973). On this basis, the domal relief of the caldera floor observed in our models may be interpreted as resulting from the forcible emplacement of experimental magma and the accomplishment of conjugated reverse faults.
Further dimensionless ratios used to analyze the modeling results are the normalized injection depth ID n and the normalized diameter of the resurgent area D n . In particular, ID n relates the position of the renewed magma chamber to the total overburden thickness (OT) above the collapsed magma chamber (ID n = ID/ OT). D n is instead obtained dividing the average diameter of the resurgent area (d r ) by the average diameter of the caldera (d c ) (D n = d r /d c ). ID n and D n ratios show a clear relationship with the type of resurgence, as outlined in the graph of Figure 9d.

Factors Controlling Surface Deformation of Domes
Piston-like resurgence generally produces little deformation on the caldera floor, unless it reactivates faults inherited from asymmetric collapse that may be located internally to the caldera. Localized resurgence affects instead part of the caldera floor, and the related experimental domes are affected by extensional faults and fractures, which in some cases may be pervasive. In these settings, extensional strain in the outer carapace of domes created by magmatic intrusion is often accommodated by outer-arc extensional faulting and fracturing (e.g., Hansen & Cartwright, 2006;Montanari et al., 2017b;Urbani et al., 2020).
In our models, the resurgent domes are deformed by normal faults showing different patterns and characteristics depending on dome dimensions and ultimately on the depth of intrusion. It is thus likely that outer-arc extension contributed to model deformation at some extent. As noted earlier, a structural pattern made up of well-defined sub-orthogonal fault segments develops for an ID of 4.5 cm (Figure 5b). The apical part of the dome becomes progressively more deformed by pervasive, closely spaced fractures and faults as ID decreases (Figures 5c and 5d). These fault-fracture networks define an overall radial pattern. In spite of the pervasive deformation, differential topography analysis reveals that sectors with differential elevation can be still identified over the experimental dome, which conceivably manifest the presence of collapsed compartments.
An analysis of fracture density was carried out on models of Series 1 showing intra-caldera doming (models LR-02-, LR-03, and LR-04). The fault pattern of these models was analyzed in great detail by mapping fault segments and each small fracture (Figure 10). We define a fracture density parameter, ρ FRAC , which is the BONINI ET AL.
10.1029/2020JB020438 20 of 32 ratio between the number of fault segments and the area of deformation induced by the resurgent dome (RIDA), and we contrast it with the applied ID ( Figure 10). The results of this analysis reveal a clear correlation between these parameters, with fracture density increasing with decreasing ID.
Such ID-dependent deformation features (i.e., fault-fracture networks and dome compartmentalization) evidently develop to accommodate the vertical and lateral growth of experimental domes. Whereas resurgent domes grow vertically ( Figure 6) in response to intrusion-driven forced folding, lateral expansion may be related to their gravitational instability and collapse. In addition, the overall centrifugal velocity gradient produced by the intrusion-bounding reverse faults should result in the stretching of the apical region of the experimental dome (see DPIV analysis of velocity vector distribution in Figure 8). Thus, this process would reasonably contribute to dome collapse, with formation of sectors with differential vertical displacement, as well as to the radial growth of resurgent domes. Deformation of experimental domes may thus be accomplished by a complex interplay among various components, particularly: (i) outer-arc extension of the growing dome, (ii) gravitational collapse of the dome, and (iii) outwards movement of inward-dipping reverse faults bounding the resurgent dome.

Role of Pre-existing Discontinuities in Caldera Resurgence
The occurrence of pre-existing discontinuities in the system has affected the intrusion process of the analogue magma at some extent. DPIV analysis applied to models with existing discontinuities show, in fact, a clear control of these features on the modalities of caldera resurgence. As noted earlier, Model DR-01 is characterized by resurgence-related radial velocity field, in which horizontal velocities are highest along the pre-existing faults and on the opposite side ( Figure 8). This behavior demonstrates that pre-existing faults represent efficient mechanical features accommodating piston-like caldera resurgence.
The role of pre-existing features can be also inferred from comparison between Model LR-02 and Model DR-02. Both models had a same ID but the final deformation pattern were rather different, with the former being characterized by a saucer-shaped sill and the latter by a dike-like intrusion (cf. Figures 5b and 6d). This behavior may be partly related to the much less quantity of experimental magma that was injected into BONINI ET AL. Model DR-02 (Table 2), because it was stopped soon after the fluid extruded at surface. However, this raises the question of which conditions did produce such a fast ascent of the magma-simulating fluid up-through the brittle overburden. The only difference with Model LR-02 was the presence of a central vertical discontinuity, and thus the reason of such a discrepancy can be sought in the rapid reactivation and exploitation of this discontinuity by the pressurized analogue magma.

Kinematics of Faulting during Asymmetric Caldera Collapse and Resurgence
Asymmetric collapse of caldera systems, due to incomplete development of ring faults or simply to differential subsidence, generates trapdoor structures (e.g., Cole et al., 2005;Lipman, 1997). Our models replicated this process and focused on the role that this geometry may play during caldera resurgence stages. As noted earlier in Section 5.2, cross sections of models corroborate the interpretation that the external normal faults reversed their movement to accommodate the lifting of the collapsed caldera block. However, most of vertical uplift was accommodated by a normal fault system delimiting the apical, and more uplifted part of the experimental trapdoor block. Such a main fault system shows a clear along-depth steepening, and develops close to the caldera rim as fault antithetic to the ring normal faults (Figures 6b and 7).
Notably, this normal fault system develops exclusively in asymmetric collapse systems (cf. Figures 11a and 11b). Interestingly, such faults dip oppositely to the block tilting and reveal their role as "lifter" of the trapdoor block, which is hinged on the opposite side of the experimental caldera (Figures 11a). The development of this fault system has also the consequence to compartmentalize the experimental caldera floor into the lifting trapdoor block and the more stable part of the caldera floor.
Finally, it might be worth noting that this normal fault system accommodated both the trapdoor caldera collapse and the subsequent resurgence of the tilted roof block without changing its kinematics. In other terms, caldera resurgence uplifted dominantly the footwall of this fault system. Model outcomes thus emphasize the role of these faults in accommodating the rotation of caldera blocks experiencing asymmetric collapse (Figures 11a).

Comparison With the Los Potreros Caldera
The LPC fault system exhibits a complex pattern with fault segments departing at high angle from the main ∼ NNW-trending Maxtaloya and Los Humeros fault segments. According to some studies, such an intra-caldera fault pattern is related to fault-bounded resurgent blocks (Norini et al., 2019(Norini et al., , 2015. The LPC fault pattern shows a remarkable similarity with Model LR-02, which is characterized by a main system of faults crosscutting the whole resurgent area and delimiting sub-orthogonal secondary faults (cf. Figures 12a  and 12b). Furthermore, the main normal fault system has a geometry and position in the experimental caldera that are similar to the WSW-dipping Maxtaloya and Los Humeros fault segments. Specifically, in both Model LR-02 and LPC, such faults are located close to the caldera rim and dip toward it. In addition, both the LPC and model LR-02 show similar fault distribution patterns, which are characterized by sub-orthogonal main frequency peaks (cf. Figures 12a and 12b). The model cross section provides further hints into the geometrical setting of this fault system. As noted earlier in Section 5.1, intrusion of the analogue magma is accommodated by steep, inward-dipping reverse faults lifting a portion of the experimental caldera block, and the model normal fault dies out at depth above a reverse fault (Figure 12c). Such geometrical relationships are comparable to those inferred by Norini   with respect to a subsurface resurgence-related reverse fault revealed by focal mechanism solutions and Formation Micro Imaging logs recorded in a geothermal well (Figure 12d). Another possibility suggested by the modeling is that the Maxtaloya-Los Humeros fault segments reactivated early normal faults associated with asymmetric (trapdoor) caldera collapse, in a similar fashion of models PR-02-PR-05 of Series 2 (see Section 5.5.4; Figures 6b, 7, and 11a). As for the presence of a central discontinuity in the model caldera (tested by Model DR-02), it did not result in a surface fault pattern similar to the natural prototype, although the discontinuity was reactivated at depth and exploited by the intruded analogue magma (see Sections 5.3 and 5.5.3). However, we cannot rule out the possibility that an eccentric rather than central intrusion point could lead to the complete reactivation of the discontinuity.
If similarity conditions are satisfied, the good correspondence between faults in model and nature may indicate a similarity in dynamic processes. Extrapolating the modeling results to nature, we speculate that the sub-angular intra-caldera fault pattern is consistent with a relatively shallow magma intrusion (scaled ID of ∼4.5 km), and not with piston-like resurgence created by renewed magma injection into the deeper chamber that originated the caldera collapse (i.e., the Zaragoza magma chamber). In addition, the eruption of Holocene trachytes, trachyandesites, and basaltic andesites around the exposed or inferred rim of the LPC (Carrasco-Núñez et al., 2017a; Figure 2) may reveal a partial reactivation of peripheral caldera faults during resurgence. This observation points to the absence of a rheological boundary at the level of the Zaragoza chamber (see Lucci et al., 2020), which has been envisaged as a heterogeneous reservoir arranged as semi-connected high-melt lenses within a partially consolidated crystal mush (Carrasco-Núñez et al., 2012).
The comparison of the LPC with experimental results can be further tested by plotting its D n and ID n ratios (see Section 5.5.1) on the graph of experimental data point distribution of Figure 9d. Unfortunately, the available data allow to calculate with confidence only D n , while the parameter ID n can be only evaluated through the correlation with the experimental best-fit curve (red dashed line in Figure 9d). Notably, the D n ratio of LPC is very similar to that of model LR2, which has been correlated to LPC on the basis of the similar resurgence-related surface fault pattern. In addition, if one considers the intersection of D n with the best-fit curve, the resulting ID n (∼0.8) indicates an ID shallower than the magma chamber that produced the caldera. This ID n value is again similar to that of model LR2 (ID n = 0.75), and thus it basically supports the initial correlation between LPC and model LR2.
The evolutionary model outlined above for the LPC integrates with that of Urbani et al. (2020), in that the ascent of shallow (<1 km depth) multiple magma bodies (emplacing as lava domes and cryptodomes) could follow structural pathways associated with caldera resurgence. This possibility may be consistent with the observation that lava domes and cryptodomes are not randomly distributed, but appear to cluster along the main intra-caldera fault zones ( Figure 2). The emplacement of these magma bodies might have also increased locally the vertical fault separation, as observed for instance at the Los Humeros fault. A corollary to the above observations is that fault-controlled shallow emplacement of magma has obvious implications for targeting geothermal fluids, which are expected to impregnate fault damage zones (e.g., Jentsch et al., 2020;Montanari et al., 2017b).
In summary, the good correlation between laboratory models and the natural fault pattern suggests a similarity in dynamic processes, and thus we propose that LPC evolution has been controlled by shallow magma resurgence.
BONINI ET AL.

Comparison With Other Resurgent Calderas and Open Challenges
The results of the current modeling revealed that caldera resurgence is generally accompanied by specific features, particularly (i) the (partial or total) reactivation of peripheral caldera faults as reverse faults, (ii) the development of conjugated reverse faults emanated from the margins of sill-like intrusions that fold the roof block in combination with the injection of analogue magma, and (iii) the development of a typical fault pattern over the surface resurgent dome that closely depends on the ID. In addition, the modeling results suggest that resurgence is not necessarily related to the recharging of the magma reservoir that produced caldera collapse but may be a more complex process associated with post caldera multistorey intrusions at various and shallower levels. Such intrusions deform smaller and localized areas of the caldera floor depending or their depths (see Section 5.5.1). Evidence for resurgence associated with polyphased and multilevel intrusions are found at Campi Flegrei caldera, where historical pulsed uplifts and bradyseisms have been modeled by repeated injection of magmas at depths as shallow as 3-4 km (e.g., Amoruso et al., 2008;2017).
In natural resurgent calderas, resurgence-related uplift is often accommodated by normal faulting along peripheral caldera faults. For instance, the analysis of seismicity that affected the Long Valley caldera, eastern California, in 1989 revealed a circular ring-normal fault that was slipping to accommodate uplift of the center of the caldera (Prejean et al., 2003). However, reverse movement of the peripheral caldera faults has been also hypothesized for some resurgent calderas, for example the Acoculco caldera complex in Mexico (Avellán et al., 2020;Figure 12e), and possibly the Xiangshan caldera in China (Jiang et al., 2005).
Reverse faults bounding intra-caldera resurgent domes have been also identified in a number of calderas (in addition to the LPC), such as the Timber Mountain caldera, Nevada (Byers et al., 1976), the Latera caldera, Italy (Tibaldi & Vezzoli, 1998), the Ischia and Pantelleria calderas in Italy (Orsi et al., 1991), the Neapolitan Yellow Tuff caldera, Italy (Corradino et al., 2021), the Long Valley caldera in California (Acocella et al., 2001), and the Reforma caldera in Baja California Sur (Pellicioli et al., 2019). Other calderas showing similar short-wavelength arching of the caldera floor might be related to similar processes (e.g., the Valles caldera in New Mexico, the Creede caldera in Colorado, and the Toba caldera in Sumatra).
The Acoculco caldera complex (ACC; López-Hernández et al., 2009) provides further hints on the relationship between ID and surface fault pattern. This feature lies some 85 km west-northwest to the LPC and shows a 18 × 16 km caldera depression controlled by pre-existing fault systems (Avellán et al., 2020). According to Avellán et al. (2020), resurgence of the caldera floor began in Quaternary (0.9-0.016 Ma), during which peripheral caldera faults were reactivated and new normal faults formed over a sub-circular resurgent area displaying a sub-orthogonal, radial pattern (Figures 12e and 12f). Such intra-caldera faults define a main ∼ E-W-trending apical graben that is interrupted in the central part of the dome (Avellán et al., 2020;Figure 12e). N-S to NNW-SSE-striking fault segments develop to the north and south with respect to the apical part of the dome, and other faults with different orientation affect more marginal sectors of the ACC (Figure 12f). The overall intra-caldera fault pattern is thus characterized by two sub-orthogonal fault systems, in a similar fashion to the LPC, even though faults as important as the Maxtaloya-Los Humeros faults are apparently missing and the portion of caldera floor affected by resurgence is seemingly larger in the ACC.
On this basis, the ACC deformation pattern may be qualitatively correlated again with Model LR-02, which shows the reactivation of peripheral ring faults and the development of apical normal faults showing a characteristic sub-orthogonal pattern (Figures 12e and 12f). Crude application of model results to the natural prototype would thus suggest a scaled ID of 4.5 km. On the basis of magnetic anomalies, Avellán et al. (2020) hypothesize that caldera resurgence was driven by renewed magma overpressure at the base of the Jurassic-Cretaceous limestone (Figure 12e), which lies approximately between 4 and 5 km depth (Calcagno et al., 2018). If this correlation is correct, then the natural ID would scale very well with that applied to Model LR-02. The good correspondence between the ID predicted from the experimental intra-caldera surface fault pattern and that inferred for the natural prototype corroborates the validity of this modeling approach (cf. Figures 12b and 12f). This correlation can be further tested using the D n and ID n ratios, as done for the LPC (section 6.1). The ACC shows a D n value slightly higher than that of LPC (and of model LR2), and the correlation with the experimental best-fit curve yields ID n ∼0.9 (Figure 9d). The latter value indicates a renewed magma ID near to the top of the early magmatic chamber, which may be consistent with the piston-like deformation proposed by Avellan et al. (2020) for the ACC. Overall, the experimental results show a satisfactory correlation with the parameters and deformation features of the investigated natural calderas (i.e., LPC and ACC).
As a final remark, natural resurgent domes often show an elongated shape in map-view (e.g., the Valles, Long Valley, Siwi, Creede, Timber Mountain, and Toba calderas; Brothelande et al., 2016;Christiansen et al., 1977;De Silva et al., 2015;Hildreth, 2017;Kennedy et al., 2012;Lipman, 2006;Phillips et al., 2007), which differ from the features modeled in this study. Smith and Bailey (1968) perceived the importance of this factor, and qualitatively compared the fault patterns of resurgent domes with that produced in an elongated dome of a clay model. Brothelande and Merle (2015) addressed experimentally this issue and investigated successfully the relationships between ID and fault pattern over resurgent domes with different elongation. However, existing caldera structures may also play some role in this process. A future challenge of analogue modeling will thus be the study of the complex interrelations among surface fault pattern, ID and caldera-forming structures in elongated resurgent domes.

Concluding Remarks
Analogue modeling has tested different conditions that may control the structural style and evolution of intra-caldera resurgence. The laboratory experiments reported here suggest the following main conclusions.
1. Models show distinct styles of deformation depending on the ID of the analogue magma (Series 1 models). Models exhibit a progression from piston-like resurgence of the whole caldera floor, for deep ID, to resurgent domes affecting an area of the caldera that decreases with decreasing IDs. Fracture density over the model resurgent dome increases as depth of intrusion decreases 2. Asymmetric caldera collapse generates the tilting of the roof block above the evacuated analogue magma reservoir, in a similar fashion to the trapdoor collapse observed in natural calderas (Series 2 models). Tilting of trapdoor block during caldera collapse is accommodated by a typical system of normal faults that bound the uplifted margin of this block and are antithetic to the nearby caldera rim normal faults 3. The presence of pre-existing faults in the pre-volcanic substratum (Series 3 models) may influence the deformation pattern of piston-like resurgence at a large extent. However, deformation outcomes produced by shallow IDs are basically similar to those of other models governed by similar conditions 4. Model results suggest that intra-caldera resurgent domes result from the forcible emplacement of analogue magma and the associated development of conjugated reverse faults lifting the interposed portion of the caldera block. Deformation of the growing resurgent dome is likely resulting from a complex in-BONINI ET AL.  (2015).

Table 3
Analogue Modeling Parameters and Scaling terplay among different parameters, including outer-arc extension and gravitational outward movement of the experimental dome, and the outwards displacement produced by the reverse faults bounding the resurgent dome 5. The experimental results appear to apply to some resurgent calderas. In particular, the fault pattern of models deformed by shallow intrusions is strikingly comparable to that of some natural volcanic complexes experiencing intra-caldera resurgence, namely the Los Potreros and partly the Acoculco caldera complexes in the eastern TMVB. Both model and nature show, in fact, a typical sub-orthogonal fault pattern, which is interpreted as evidence for similarity in dynamic process. Application of model results would suggest that surface deformation of caldera floor may be induced by renewed magmatic pressure produced at ∼4.5 km depth

Appendix A: Model Scaling
Definition of geometrical, kinematical, rheological, and dynamical model to nature ratios is essential for evaluating appropriately the modeling results and to compare them with the natural prototype (Hubbert, 1937;Ramberg, 1981;Weijermars & Schmeling, 1986). Models are conveniently scaled such that 1 cm in the model corresponds to 1 km in nature, involving a model to nature geometrical ratio l* = l m /l n = 10 −5 (subscripts m and n refer to models and nature, respectively; see Table 3). Thus models simulate a brittle overburden that scales to a thickness of 6 km, and a caldera collapsed area that has dimensions similar to the LPC. The model to nature ratio of stress is σ* = σ m /σ n = ρ*g*l* = 5.4 × 10 −6 , considering a brittle sedimentary cover with average density ρ n = 2,600 kg m −3 (ρ* = ρ m /ρ n = 0.54), and g* = 1 (models have been deformed under a normal gravity field). A similar density ratio ρ* can be obtained for granitic to andesitic melts with density varying approximately between 2,200 and 2,500 kg m −3 (Murase & McBirney, 1973). σ* is very sensitive to the chosen length ratio l*, and the assumption of a Mohr-Coulomb failure criterion for shallow crustal rocks implies that cohesion c, having the dimensions of stress, should be scaled with the same ratio, particularly c m = c n σ*. Mean values of cohesion of shallow crustal rocks span over a range of 5-40 MPa depending on the rock type (Hoshino et al., 1972;Jaeger & Cook, 1969). However, in magmatic systems rocks are intensely fractured and deeply altered by the circulation of hydrothermal fluids, whereby cohesion may be significantly lowered (Merle, 2015, and references therein). For an average cohesion of 10 MPa, the correspondent analogue material should have a cohesion of 54 Pa, which is again of the same order of magnitude of our sand mixture (∼10 Pa). In nature, sedimentary rocks have a coefficient of friction μ n ranging from 0.6 to 0.85 (Byerlee, 1978), and thus our experimental overburden falls in the upper limit of this range (μ m = 0.83).
Dynamic similarity requires that model to nature ratios of different forces (e.g., gravitational, inertial, viscous, elastic, and frictional) are constant (Hubbert, 1937;Ramberg, 1981). This condition implies that the ratio between any two different forces must be identical in both model and natural prototype (Hubbert, 1937;Ramberg, 1981). The evaluation of dynamic similarity is thus based on this condition, and can be tested through the comparison of dimensionless numbers given by the ratio of characteristic forces acting on models and prototype.
One way to determine dynamic similarity is to assume that models are scaled to nature, and then compute the mean bulk viscosity of natural magmas and compare it with actual values (e.g., Bonini, 2003;Sokoutis et al., 2000). The Ramberg number (R m , Weijermars & Schmeling, 1986) is optimal for this purpose, as this dimensionless number represents the ratio of gravitational force to viscous force acting on a viscous layer, which is given as Ramberg (1981): where g is the gravitational acceleration, ε is the strain rate, and ρ, η, and l are density, viscosity and characteristic length (thickness and width) of the viscous layer.
We focus our analysis on model resurgence, taking the simpler case in which the analogue magma is re-injected into an existing magma chamber (piston-like resurgence). The experimental Ramberg number exp m R can be calculated for model resurgence considering ρ m = 1,190 kg m −3 , g m = 9.81 m s −2 , l m = 0.12 m (i.e., the diameter of the analogue magma chamber), η m = 17 Pa s, and ε m = 1.85 × 10 −3 s −1 . The strain rate ε m has been evaluated for a Newtonian fluid undergoing near pure shear deformation, and approximated as    / x v x (Ramberg, 1980), where v is the applied injection velocity (80 cm h −1 = 2.2 × 10 −4 m s −1 ) and x is the width of the deformed zone in the model (i.e., the magma chamber width, 0.12 m). The above parameters yield exp m R = 4.4 × 10 4 . In case of injection velocity of 40 cm h −1 (1.1 × 10 −4 m s −1 ; Table 2) exp m R = 9.03 × 10 4 .
A mean bulk viscosity  n is estimated for the natural magma assuming that the experimental Ramberg number exp m R equals that of the natural prototype: n n n n m n g l R (2) Dynamic scaling of magmatic systems is difficult because many key parameters are often unknown, particularly, density, viscosity, and ascent velocity of magma (Merle, 2015). We therefore assume average values, specifically ρ n = 2,400 kg m −3 (an average density of andesitic magma at emplacement conditions; Murase & McBirney, 1973), g n = 9.81 m s −2 , l n = 12 × 10 3 m (the scaled magma width), and ε n = 1.67 × 10 −8 s −1 , assuming ε n = v n /l n and an average velocity of magma emplacement v n = 2 × 10 −4 m s −1 (e.g., Merle, 2015). The viscosity (η n ) estimated from eq. (2) for the natural melts is of the order of 3.8 × 10 11 Pa s for the highest injection velocity (80 cm h −1 ) and 1.8 10 11 Pa s for the lowest injection velocity (40 cm h −1 ), respectively. These viscosities fall approximately in the lower range for acid magmas, which may vary from 10 11 -10 12 Pa s up to 10 15 Pa s for crystal-rich magmas (Merle, 2015, and references therein). The obtained magma viscosity is thus enough similar to actual natural magma viscosities to validate our initial assumption of dynamic similarity between models and natural prototype. Since both natural and analogue overburden share similar parameters of brittle behavior (see Table 3), we can extrapolate our conclusion of dynamic similarity to the whole caldera system.
The ratio of ascent velocity of magma is given by the relationship v* = ε* l* (Merle & Vendeville, 1995), which can be used to obtain the ratio of time t*. Since strain rate and time ratios are related as ε* = 1/t*, the ratio of time can be written as t* = l*/v*, and thus t n = (t m v*)/l*. Recalling that in these models it has been assumed v* = 1.1 for v m = 2.2 × 10 −4 m s −1 (v* = 0.55 for v m = 1.2 × 10 −4 m s −1 ) and l* = 10 −5 , it implies that one hour in the model corresponds to 12.6 years in nature (6.3 years for v m = 1.2 × 10 −4 m s −1 ). Finally, the combination between the assumed natural ascent velocity (2 × 10 −4 m s −1 ) and the scaled diameter of the injection pipe (400 m) yields a scaled volumetric magma flux of 0.79 km 3 yr −1 .

Data Availability Statement
All the data supporting this research are available in the text, in the Appendix and in the in-text citation Maestrelli et al. (2020) free to download at the following link http://doi.org/10.5281/zenodo.3882130 under Creative Commons Attribution 4.0 International Licence.