Front. Earth Sci. Frontiers in Earth Science Front. Earth Sci. 2296-6463 Frontiers Media S.A. 1535083 10.3389/feart.2025.1535083 Earth Science Original Research Stealthy magma system behavior at Veniaminof Volcano, Alaska Li et al. 10.3389/feart.2025.1535083 Li Yuyu 1 * Gregg Patricia M. 1 Lu Zhong 2 Wang Jiahui 2 1 Department of Earth Science and Environmental Change, University of Illinois Urbana Champaign, Urbana, IL, United States 2 Department of Earth Sciences, Southern Methodist University, Dallas, TX, United States

Edited by: Gilda Maria Currenti, National Institute of Geophysics and Volcanology (INGV), Italy

Reviewed by: Vyacheslav M. Zobin, University of Colima, Mexico

Yan Zhan, The Chinese University of Hong Kong, China

*Correspondence: Yuyu Li, yuyuli3@illinois.edu
10 06 2025 2025 13 1535083 26 11 2024 07 04 2025 Copyright © 2025 Li, Gregg, Lu and Wang. 2025 Li, Gregg, Lu and Wang

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Although Veniaminof Volcano in Alaska experiences frequent eruptions and has eight permanent seismic stations, only two of the past 13 eruptions have had precursory signals that prompted a pre-eruption warning from the Alaska Volcano Observatory (AVO) since 1993. Seismic data from Venianimof indicate that most eruptions from 2000 to 2018 do not coincide with increased seismicity. Additionally, analyses of InSAR data available from 2015 to 2018 which covers the pre-, syn-, and post-eruption periods of the 2018 eruption do not show clear signs of deformation. The systemic lack of systematic precursory signals raises critical questions about why some volcanoes do not exhibit clear unrest prior to eruption. Volcanoes that erupt frequently without precursory signals are often classified as “open” systems with magma migrating through an open network to eruption, rather than pausing at a shallow reservoir. However, the precursory signals, or lack thereof, from a small or deep closed magma system may be difficult to observe, resulting in a stealthy eruption mimicking the behavior of an open system. In this study, we utilize finite element, fluid injection models to investigate a hypothetical closed magma system at Veniaminof and evaluate its ability to erupt with no observable early-warning signals. Specifically, a series of numerical experiments are conducted to determine what model configurations lead to stealthy eruptions – i.e., producing ground deformation below the detection threshold for InSAR (<10 mm) and developing no seismicity, yet resulting in tensile failure which will promote diking and eruption. Model results indicate that the primary control on whether eruption precursors from deformation and seismicity will be present are the rheology of the host rock and the magma flux, followed by the secondary control of the size of the magma chamber, and then its depth and shape. Volcanoes with long-lived thermally mature magma systems with moderate to small magma reservoirs are the most likely to exhibit stealthy behavior, with the smallest systems most likely to fail without producing a deformation signal. This result is likely because small, deep magma systems produce minimal surface deformation and seismicity. For stealthy volcanoes like Veniaminof and others in Alaska (e.g., Cleveland, Shishaldin, Pavlof) and around the world, understanding the underlying magma system dynamics and their potential open vs. closed nature through numerical modeling is critical for providing robust forecasts of future eruptive activity.

Aleutian volcanoes Veniaminof eruption volcano forcasting magma system finite element model rheology section-at-acceptance Volcanology

香京julia种子在线播放

    1. <form id=HxFbUHhlv><nobr id=HxFbUHhlv></nobr></form>
      <address id=HxFbUHhlv><nobr id=HxFbUHhlv><nobr id=HxFbUHhlv></nobr></nobr></address>

      1 Introduction

      Volcano monitoring and forecasts are crucial to evaluate and mitigate the potential risks and socioeconomic impacts of volcanic activity. Seismicity, ground deformation, and gas emissions are commonly used as precursory signals and play a critical role in providing early warning of volcanic activity and impending eruption (e.g., Segall, 2013; Voight et al., 1998). However, sometimes observations of deformation and/or increase in seismic activity do not correlate with the timing of eruptions (Biggs et al., 2014; Roman and Cashman, 2018). Evaluating volcanoes that lack eruption forecasting signals is necessary to understand the characteristics of stealthy magma systems and develop strategies to mitigate their impacts. Veniaminof, an active caldera system on the central Aleutian Arc of Alaska, exhibits the characteristics of these stealthy eruptions (Figure 1). Although Veniaminof experiences frequent eruptions and currently has eight permanent seismic stations, only two of the past 13 eruptions since 1993 have had precursory phases that prompted a pre-eruption warning from the Alaska Volcano Observatory (AVO; Cameron et al., 2018). Previous research has indicated several hypotheses for why a volcano may erupt without warning, including: magma reservoirs that are too deep to create detectable inflation signals during magma recharge and the migration of eruptive products from a deep mid-crustal reservoir that only transit through the shallow magma storage systems (Grapenthin et al., 2013); magma flux that is below the detectible level; and open volcanic systems that have open vents and frequent degassing resulting in a lack of detectable pre-eruptive ground uplift precursors (Chaussard et al., 2013; Biggs et al., 2014). However, the spectrum between open-vent and closed-system unrest is likely not dichotomous but rather relative, with volcanoes experiencing variations between open and closed modes of unrest (Rose et al., 2013; Acocella et al., 2023). Well-recognized open-vent volcanoes such as Villarrica (Chile), Stromboli (Italy), and Galeras (Colombia) may represent archetypal cases. In contrast, Veniaminof may occupy a position that incorporates elements of both models. The precursory signals from a small or deep closed magma system may be difficult or impossible to distinguish from open vent unrest, as both result in eruptions without observable precursors.

      Veniaminof Volcano located in the Aleutian Arc, Alaska. (A) Coastline map of the volcanoes of Aleutian Arc. Red box represents the location of (B). (B) Topographic map of Mt. Veniaminof showing its summit (red triangle) and the locations of eight permanent seismic stations maintained by AVO (yellow squares).

      The goal of this study is to investigate the magmatic processes in closed systems that lead to stealthy eruptive behavior characterized by the absence of seismicity (i.e., shear failure) and observable ground deformation, while still resulting in tensile failure leading to dike initiation and eruption. In other words, we aim to investigate how a closed magma system may behave in a way that mimics open system unrest. We use seismicity and ground deformation observations from Veniaminof in conjunction with finite element method (FEM) modeling (Le Mével et al., 2016; Gregg et al., 2012) to address specific questions regarding volcanic activity that does not produce precursory signals of an imminent eruption: 1) What magma system parameters have the greatest effect on observable eruption precursory signals? 2) What parameters are necessary to trigger stealthy eruptions from closed magma systems? 3) What do these findings indicate about the magma system at Veniaminof–such as the associated inter-eruptive behavior, the stability of the system, and its potential for eruption? Axisymmetric, 2D fluid injection models are utilized to identify the configurations that produce ground deformation signals below the detection threshold of InSAR (<10 mm) and no seismicity, while requiring a minimum magma flux needed to produce tensile failure and drive eruption. Elastic models with non-temperature-dependent and temperature-dependent rheology are applied to Veniaminof Volcano to evaluate the overall effect on observation of seismicity and ground deformation. To constrain the cases that result in stealthy eruptions, ones that occur without visible precursory signals, our particular focus is the effects of parameters such as the size, geometry, depth of the magma chamber, and the flux of magma injecting into the conduit and chamber, on the displacement, tensile rupture, and Mohr-Coulomb failure related to the observable precursor signals of InSAR deformation and seismicity.

      2 Geologic background

      Mt. Veniaminof is an ice-clad, basalt-to-dacite stratovolcano located on the Alaska Peninsula, ∼750 km southwest of Anchorage (Figure 1). It is one of the largest caldera-forming volcanoes of the Aleutian arc and is comprised of a 10 km diameter caldera surrounded by 600 m wall scarps formed during caldera-forming eruptions 3.7 kya with an estimated volume of ∼350 km3 (Miller et al., 1998; Bacon et al., 2007). Veniaminof is one of the most active volcanoes in the Aleutians and has erupted at least 15 times in the past 200 years from the intra-caldera cone (∼300 m in height). Eruption compositions are primarily basaltic to basaltic andesite, resulting in the formation of minor tephra deposits and small lava volume effusive flows mostly within the boundaries of the caldera (Loewen et al., 2021). Veniaminof is surrounded by extensive pyroclastic-flow deposits from two Holocene caldera-forming eruptions, nine kya and 3.7 kya, and lahar deposits (Miller and Smith, 1987; Bacon et al., 2007). Historical eruptions have produced small-volume basaltic to andesitic deposits from an intra-caldera cone, while the post-3.7 kya dacitic pumice-fall deposit represents explosive silicic eruptions (Fournier and Freymueller, 2008; Bacon et al., 2007). Previous work indicates that the subsurface magma system of Veniaminof consists of a shallow reservoir comprised of a region of segregated felsic magma on top of a crystal mush column where basaltic magmas occasionally intruded (Bacon et al., 2007). The observation of anhydrous phenocryst assemblages and whole-rock geochemistry indicates a dry, reduced, and high-temperature magmatic system (Miller and Smith, 1987).

      Since the 1830–1840s, 19 eruptions have been documented (Waythomas et al., 2022), with most occurring without clear early warning, such as the 2004, 2005, 2006, and 2008 eruptions, which were classified as “detect only” by the Alaska Volcano Observatory (Cameron et al., 2018). The most recent eruptions in 2021 and 2018 also lacked distinct precursory signals. The 2021 eruption (VEI 1) was not detected until 3 days after its onset, when sulfur dioxide emissions were observed, followed by increased surface temperatures, an ash cloud, and an ash plume from Cone A, which consisted of small explosions and minimal lava effusion within the glacial melt pit 1 km east of Cone A (Waythomas, 2021; Smithsonian Institution Global Volcanism Program, 2021; Orr et al., 2024). While the 2021 eruption has not been analyzed in detail, recent investigations of the 2018 eruption provide insights into the magma system. In fall of 2018, a relatively long-duration eruptive event occurred (VEI 2), lasting from September 4th to December 27th 2018. The strombolian eruption generated a lava flow 800 m long, and minor ash emissions (Smithsonian Institution Global Volcanism Program, 2021). From September to December, the explosivity of the eruption increased producing a continuous ash plume of 400 km on November 21st (Loewen et al., 2021). In addition to the observation from the local seismic networks and satellite images, tephra and lava flow samples have been collected (Wallace et al., 2020; Loewen et al., 2021). Preliminary analyses, including glass geochemistry and whole-rock studies, indicate basaltic andesite products initiating from ∼16 km depth with consistent compositions throughout the eruption period (Loewen et al., 2021). While 2018 eruption exhibited low-level volcanic tremor, slight surface temperature increases, and minor steam emissions immediately before its onset, these signals resembled past non-eruptive unrest and occurred too close to onset for reliable warnfing (Waythomas et al., 2022).

      Changes in seismicity rate can be a precursory signal of imminent eruption. However, Veniaminof eruptions during 2000–2021 either mismatch changes in cumulative seismicity or have no concurrence with rapid increases in seismicity (Figure 2; Power et al., 2019). Additionally, several increases in seismicity occurred between 2008 and 2013 and no eruptions occurred during this 5-year period. Changes in the cumulative seismicity may correspond with eruptions in 2004, 2008, and 2015, but there is a dearth of shallow earthquakes. Overall, the number of the shallow earthquakes associated with volcanic unrest/eruption at Veniaminof is minimal. These observations indicate that seismicity is not a reliable precursory signal for Veniaminof eruptions, leading to the question of whether Veniaminof is an open system (e.g., Chaussard et al., 2013) or a closed, stealthy magmatic system. To determine the dynamics underlying Veniaminof’s activity additional characteristics of its magma reservoir must be examined.

      Earthquake hypocenters from 2000 to 2021 near Veniaminof Volcano, Alaska. (A) Earthquake epicenters within 20 km of the summit color coded by magnitude. (B) Longitude versus hypcentral depth. The red triangle indicates the location of the Veniaminof summit. (C) Hypocentral depth versus time and cumulative number of located earthquakes (red line). Green regions indicate eruptive periods. Most eruptions between 2000 and 2021 are not proceeded by increased earthquake activity, while no eruptions between 2008 and 2021 occurred when several increases in earthquake activity were detected. (Power et al., 2019 and personal communication from Power, J.).

      3 Methods 3.1 Persistent Scatter SAR interferometry (PS-InSAR)

      To observe the spatiotemporal evolution of the ground deformation before and throughout the 2018 Veniaminof eruption (September to December 2018), Persistent Scatter SAR Interferometry (PS-InSAR) (Ferretti et al., 2000) was applied on the Sentinel-1 data from 2015 to 2018. The summer data (June to October) are selected for each year to obtain the best coherence between images, as the images are impacted by the ice and snow coverage at the caldera of Veniaminof. The SAR images are C-band (wavelength = ∼5.56 cm) Single-Look-Complex (SLC) Level 1 products of Sentinel-1A from ascending tracks, using IW swath mode with VV polarization and incidence angle 20°–46°.

      StaMPS was used to perform time-series analysis of the SAR acquisitions (Hooper et al., 2004; Hooper, 2006; Hooper et al., 2007; Hooper et al., 2007; Hooper, 2008). The topographic phase was removed from the interferograms with Shuttle Radar Topography Mission (SRTM) one arc-sec (30 m) data. Orbit error and atmospheric phases were estimated and subtracted from differential interferograms during the StaMPS processing.

      3.2 Finite element methods

      To investigate host rock stability in response to the ground deformation caused by inflation of magma chamber, we use a thermomechanical Finite Element Method (FEM) modeling approach. We build upon previous numerical experiments (Gregg et al., 2012; Grosfils et al., 2015; Le Mével et al., 2016; Gregg et al., 2018) utilizing COMSOL Multiphysics 6.0 to model an inflating magma reservoir at Veniaminof Volcano (Figure 3). The COMSOL modules utilized in this investigation include Heat Transfer Module, Fluid Transfer Module, Nonlinear Materials Module, and Structure Mechanics Module. The fluid-structure interaction, time-dependent model was constructed following the setup of Le Mével et al. (2016), to solve for temperature, mass transport, and stress and strain in response to magma flux into a shallow reservoir (see Supplementary Tables 1, 2 for a list of parameters and variables). Outputs from the numerical experiments provide constraints on the failure behavior of the host rock surrounding a reservoir.

      The 40 by 16 km two-dimensional model space is constructed with a symmetrical central axis, roller boundary conditions are implemented along the side and bottom, and the ground surface is free boundary including an approximate topography of the edifice and caldera of Veniaminof (Figure 3). An existing magma reservoir is assumed with a connecting conduit to accommodate magma transporting from a deeper source (>13 km). The initial magma reservoir is varied in size from 1–6, with Size one being the largest volume and Size six being the smallest, and shape (spherical, oblate, and prolate ellipsoids). The surrounding host rock is modeled as a linear elastic material. For additional details of the model implementation please see the Supplemental Methods.

      Finite element method (FEM) model set up for Veniaminof. The fluid injection approach of Le Mével et al. (2016) is modified to investigate the magma system at Veniaminof and its ability to erupt with no geophysical precursors. An axisymmetric 2D model is set up with side and bottom roller boundaries and the free boundary on surface, including the approximate topography of Veniaminof volcano. Magma intrudes upward from a deep-seated source (>13 km) via a conduit into a shallow reservoir. Volumetric growth of the magma chamber is caused by the magma injection. All the following displacement results (Figures 57) are taken as the edge of the caldera (yellow dot in this figure).

      4 Results 4.1 InSAR-derived ground deformation

      The time series of observed ground deformation from summer 2015 to 2018 (before and during 2018 eruption) covering the area of volcano edifice of Veniaminof generated by PS-InSAR are provided in Figure 4. The line of sight (LOS) displacements are in the range of ±3 cm including motion towards the satellite (uplift) and away from the satellite (subsidence). The time series deformation results at each coherent pixel were presented in Supplementary Figure A.1. The uplift and subsidence of a deforming volcano are considered to be the result of the movements of subsurface magma related to the unrest and eruption activities (Dzurisin, 2003). However, the time series of SAR images shows no clear trend. Although our data processing involves the removal of SAR images with high noise, the displacement signals are still ambiguous, and suggest the signals observed in these images are likely due to the atmospheric effects. This lack of an observable volcanic ground deformation signal may indicate that the inflation signal caused by volcanic activity is being masked by the atmospheric effects and below the detection threshold. The detection threshold depends on number of SAR images (more images generate less uncertainty), the atmosphere condictions and topography of the study area, and the algorithms of InSAR approaches (Hanssen, 2001). In the case of Veniaminof, 1–2 cm is the approximate uncertainty of displacement. That is to say, if the precursory signal of inflation prior to the 2018 eruption is below this threshold, a clear trend may not be visible in the time series. Additionally, ground deformation at Veniaminof may be masked by the glacier located in the center of Veniaminof’s caldera. In this study we seek to link observations of ground deformation based on the ongoing InSAR data with seismicity to evaluate the stealthy nature of Veniaminof to determine the nature of its magma system that is evolving below the threshold of seismic or ground deformation to be observed as eruptive precursors.

      Ground deformation and velocity from InSAR time series from 2015 to 2018. (A,B) are pre-eruption (June 2015 to August 2018) and co-eruption (September 2018 to October 2018) ground velocity maps generated by PS-InSAR. P1 to P4 are selected locations from near edge to away from edge. The time series of ground deformation of these points are plotted in (C-F). Time series of the ground displacement (C-F) for four locations at P1, P2, P3, and P4, respectively.

      4.2 Sensitivity tests in FEM-Based simulations

      The primary goal of our numerical experiments is to determine the model configurations needed to produce stealthy eruptions - in other words, magma system states that produce a ground deformation signal that remains below the detection threshold for InSAR and results in little to no seismicity (as calculated by shear failure throughout the model space) while maintaining a minimum flux of magma needed to produce the eruptive volume recorded for the 2018 eruption of Veniaminof Volcano. As with other active systems (e.g., Katla and Grímsvötn in Iceland; subglacial volcanoes in Antarctica; Mount Spurr and Westdahl in Alaska) the center portion of Veniaminof is masked from deformation observations due to summit glacier/snow. As such, modeled deformation results are taken from the caldera rim coinciding with where the InSAR observations are coherent.

      The numerical results illustrate how observable ground deformation and the stability of magma system are controlled by the underlying magma reservoir characteristics including size, shape, depth, and magma flux (Figures 57). First, as is to be expected, shallower magma chambers are associated with higher displacement (Figure 5). Second, higher displacement is directly associated with larger chamber sizes (e.g., Size 1 and Size 2). However, for the largest reservoirs, there is also a dependence on the specific depth, shape, and flux combination as Size 2 can produce slightly higher flank displacements than Size 1 for non-temperature dependent prolate and spherical models. Third, in most cases, oblate chambers are predisposed to maximal displacement, in contrast to the minimal displacements associated with prolate chambers. Note, the difference between the prolate and oblate geometries may be due to the chosen observation spot on the edge of the caldera rather than the center of the summit. An oblate shape will generate more upward displacement as its elongated horizontal shape extends towards the edge, whereas prolate geometries result in greater deformation in the center. In addition, because the deformation is calculated at the caldera edge, the displacement does not always increase with shallower magma chamber depths (Figure 5). For example, the highest displacement value for oblate reservoirs is from −4 to −6 km (Figures 5A,B). Fourth, as is to be expected, a direct correlation exists between magma flux and displacement, with the flux rate being the most crucial factor controlling ground deformation and the timing of failure (Figures 6, 7). Negative values of displacements (<0 mm) are generated when the flux rate is <0.01 m3/s, i.e., 0.0003 km3/yr. To visualize the relationship, Supplementary Figures A2, A3 provide additional model results for flux values from 0.01 to 0.1 m3/s with steps of 0.01 m3/s to mimic the case of the flux leading up to Veniaminof’s 2018 eruption.

      Magma chamber depth vs. ground deformation. Depth is to the surface of volcano edifice with summit of ∼2 km, while ground deformation is at the edge of summit caldera (see the edge location in Figure 3) caused by inflating magma chamber. Flux is taken assumed to be 0.0507 m3/s, estimated by the volume of volcanic products produced during the 2018 eruption of Veniaminof. Magma chamber size 1 to 6 (from largest to smallest) represents the scenarios where the eruptible volume of the chamber is 0.16%, 1%, 5%, 10%, 20%, 50% of the overall reservoir volume (see details in Supplementary Table 3). The yellow shaded area indicates displacement from 0 to 10 mm, the detectable threshold of Sentinel one data. The dashed line (20 mm) marks the approximate upper limit of the detectable threshold used in Figure 7. Experiments in (A-C) assume a non-temperature dependent models, while (D-F) assume a temperature dependent rheology. Y = yes, N = no, Sm = small, Shear = shear failure, and Tensile = tensile failure. “Tensile Y” indicates tensile failure occurred along the reservoir or conduit in the experiment (assuming a tensile strength of 0 MPa as a conservative minimum).

      Magma flux vs. ground deformation. Ground deformation at the edge of summit caldera caused by inflating magma chamber, with a constant depth = −6 km (depth to the surface of volcano edifice with summit of ∼2 km height). Flux = 1e−5, 5e−5, 1e−4, 5e−4, 1e−3, 5e−3, 0.01, 0.05, 0.1 m3/s. Magma chamber size and shape parameters are described in Supplementary Table 3. The yellow shaded area indicates displacement from 0 to 10 mm, the detectable threshold of Sentinel one data. The dashed line (20 mm) marks the approximate upper limit of the detectable threshold used in Figure 7. The gray shaded region indicates models results of subsidence at the flank. Experiments in (A-C) assume a non-temperature dependent models, while (D-F) assume a temperature dependent rheology. Y = yes, N = no, Sm = small, Shear = shear failure, and Tensile = tensile failure. “Tensile Y″ indicates tensile failure occurred along the reservoir or conduit in the experiment (assuming a tensile strength of 0 MPa as a conservative minimum).

      The parameter space for stealthy eruptions. Plotted results include models that produce <20 mm of surface displacement and little to no shear failure, but exhibit tensile failure surrounding the magma reservoir. (A,B) are the results for shallow magma reservoirs (1–4 km depth below the summit) that produce stealthy eruptions (shaded gray), whereas (C,D) illustrate the regions of the parameter space for deeper reservoirs (>6 km depth below the summit) that produce stealthy eruptions (shaded green). The overlap in parameter space is indicated by the black lines for ease of comparison. See Supplementary Table 3 for the specific values of the radius of magma chamber size (1–6) for different shapes. CF = Mohr-Coulomb failure, TS = tensile stress that is over the critical threshold, and disp. = displacement.

      Our study indicates that the key trends apply in both non-temperature dependent and temperature dependent rheological conditions: 1) high flux typically leads to increased displacement and a higher chance of tensile and Mohr-Coulomb failure; 2) shallow depth is associated with high displacement but has very little impact on tensile and Mohr-Coulomb failure; and 3) larger chambers usually exhibit more displacement but are less prone to cause tensile and Coulomb failures.

      5 Discussion 5.1 Factors controlling precursory observations

      Temperature-dependent models have a lower calculated displacement and, thus, a lower likelihood of producing tensile and shear failure. Therefore, to reach the same displacement and threshold for failure, temperature-dependent models require a higher magma flux rate (Figures 5, 6). Specifically, initiating tensile failure to catalyze eruption requires a higher flux rate than the non-temperature dependent models (i.e., low fluxes of 0.01 m3/s will not suffice). Figure 7 illustrates different combination of factors and their influence on stealthy eruptions.

      1) High flux and large chamber size: high flux combined with a large size tends to produce significant displacement, excluding the likelihood of a stealthy eruption, although Mohr-Coulomb failure may not necessarily be high. This is attributed to the positive association between both high flux and large size with displacement. In temperature-dependent models (Figures 7B,D), a flux of 0.10 m3/s with large chamber size remains stealthy, and observable eruptions are beyond the plotted range due to the temperature-dependent rheology shifting the upper flux limit for stealthy eruptions rightward.

      2) High flux and small chamber size: this combination typically results in high Mohr-Coulomb failure, making seismic signals observable and unambiguous. High flux and small size jointly increase the affected subsurface region of Mohr-Coulomb failure. Yet, this combination contributes to initiate tensile failure, facilitating eruption onset, and displacement may not be significantly high, as the small size counteracts the effect of high flux. This “observable eruption” scenario is illustrated in the unpopulated areas in Figures 7A,C. With temperature dependent rheology (Figures 7B,D), the upper flux limit of “stealthy eruption” shifts rightward, accommodating higher flux values without substantially increasing Mohr-Coulomb failure.

      3) Low flux and large chamber size: a low flux coupled with a large chamber size typically prevents tensile failure (i.e., reducing likelihood of triggering of eruption) and Mohr-Coulomb failure, as both low flux and large size decrease the probability of such occurrences, although this might result in low displacement due to low flux. In non-temperature-dependent models, a low flux (0.01–0.02 m3/s) combined with a large chamber size remains effective, indicating temperature-dependent wall-rock properties amplify the combined effects of low flux and large size, thereby diminishing the likelihood of tensile failure initiation.

      4) Low flux and moderate to small chamber size: in this scenario, tensile failure and Mohr-Coulomb failure are unlikely. Although displacement is below the threshold for detection, the absence of tensile failure suggests that an eruption will not occur. While small size increases the likelihood of tensile failure, the predominant influence is the low flux.

      Therefore, the primary determinants on eruption precursors from deformation and seismicity are the rheology of the warm wall rock and the magma flux, followed by secondary parameters of the size of the magma chamber, and then its depth and shape. Essentially, a long-lived system with ample thermal input to warm the rheology has a greater parameter space that will produce stealthy eruptions without precursory signals.

      Our models, both temperature dependent and non-temperature dependent, are elastic and do not account for the viscosity of the wall rock. A viscoelastic rheology typically results in greater deformation compared to purely elastic models, as the viscous component allows for more prolonged and extensive deformation under stress before reaching failure (Bonafede and Ferrari, 2009; Del Negro et al., 2009; Gregg et al., 2012; Hickey et al., 2013; Newman et al., 2001). However, an elastic model is appropriate when the loading time is shorter than the host rock’s relaxation time, and particularly appropriate for Veniaminof which has an inter-eruption time interval of ∼5 years which should experience negligible effects due to viscous relaxation (Zhan and Gregg, 2019).

      5.1.1 Long-lived vs. transient magma systems: When does temperature matter?

      Distinguishing between transient and long-lived magmatic systems is essential for understanding volcanic behavior and associated hazards. Long-lived magma chambers exhibit sustained activity over protracted periods, often with complex, multi-tiered magma chambers, and extend through the crust and comprise heterogeneously distributed melt, crystals, and exsolved volatiles (Cashman et al., 2017). In contrast, transient magmatic systems typically feature ephemeral magma chambers that undergo rapid modifications influenced by processes such as volatile degassing, which critically impact magmatic overpressure and consequently drive volcanic eruptions (Mittal and Richards, 2019). The lifespan of a magma chamber can vary significantly depending on several factors, including the tectonic setting, the composition of the magma, and the dynamics of the magmatic system (Gualda et al., 2012; Cooper et al., 2017; Turner and Costa, 2007), and eruptible portions within long-lived chambers typically last from centuries to hundreds of thousands of years (Karakas et al., 2017).

      Veniaminof displays some characteristics of a long-lived system, such as sustained activity over millennia, but the behavior of the 2018 eruption also aligns with aspects of a transient system. Veniaminof’s eruption styles vary widely, ranging from effusive to explosive with a history of sustained volcanic activity characteristic of a long-lived system; however, historical eruptions demonstrate simultaneous explosive and effusive activity from separate vents, which could also be indicative of a transient magmatic system (Waythomas, 2021). The eruptive behavior of 2018, evolving in explosivity over time (Loewen et al., 2021; Bennington et al., 2018), might also be consistent with a transient magmatic system. Moreover, the potentially small magma chamber volume of ∼0.8–5 km3 and relatively low flux rate indicated by our models, as well as the short repose time between its frequent eruptions may suggests a transient magmatic system dynamic. Veniaminof might not maintain a continuous, high level of activity due to a consistent supply of magma, but rather due to more episodic activity, where eruptions occur in response to the accumulation of sufficient magma to initiate failure.

      In long-lived magma systems with substantial heat and material influx, the thermal state of the host rock becomes a critical factor for assessing reservoir stability (Gregg et al., 2012; Jellinek and DePaolo, 2003). For example, simulations of magma injection into a long-lived thermally primed host rock indicate that models excluding temperature-dependent elastic moduli and viscosity fail to reproduce realistic stress evolution (Gregg et al., 2018) or observed deformation level (Le Mével et al., 2016). While the inclusion of temperature-dependent rheology is an important factor for estimating the stress evolution when investigating long-lived system dynamics, the non-temperature dependent models might be applicable to Veniaminof given the potentially low flux and a small chamber size estimated by the volume of its volcanic products. The influence of the host rock’s thermal state might not be as pronounced when compared to volcanoes with larger, long-lived systems. Nevertheless, including a temperature-dependent host rock in models for Veniaminof may be appropriate. Veniaminof may represent an intermediary scenario within the constraints of parameters tested by our models, embodying a magmatic system that merges aspects of both end-member cases, cold host rock rheology (non-temperature dependent models) and warm host rock rheology (temperature dependent models).

      5.2 Stealthy vs. observable magmatic system

      Veniaminof Volcano may exemplify a stealthy eruption scenario, confirming that such eruptions are indeed feasible, an eruption characterized by the absence of observable precursory signals including seismicity and ground deformation, but sufficient tensile stress to initiate an eruption without the observable warnings.

      5.2.1 Comparison of stealthy and observable closed magma systems

      Our models are constructed based on closed volcanic systems and indicate the possibility of two types of eruptions: the “observable” eruption and the “stealthy” eruption. As depicted in Figure 8A, the “observable” eruption encompasses various activity stages. During the pre-eruptive stage, the influx of new magma or the exsolution of volatiles leads to reservoir pressurization and volcanic edifice inflation. This inflation, both in distance and elevation on the surface, can be detected using methods such as InSAR or GPS (Figure 8A1). Additionally, during this stage, the wall rock experiences stress due to magma ejection from the deep source. As a result, rupture of the wall rock occurs around the chamber and conduit region, accompanied by a significant number of volcano-tectonic earthquakes corresponding to shear failure. When the pressure surpasses the crustal strength and the tensile stress on the wall rock exceeds the critical threshold of 0 MPa (chosen as an end-member, minimum value), magma can ascend through a conduit/dike toward the surface, initiating the volcano’s eruption (Figure 8; Supplementary Figure A2). Once the eruption commences, during its co-eruptive stage, the removal of magma and gas from the reservoir causes the magma chambers to contract and the pressure to decrease. Consequently, deflation surface signals can be observed, unless there is additional recharge during the eruption. Additionally, in the usual case, volcano-tectonic earthquakes disperse during this stage. During the post-eruptive (cooling) stage, characterized by the deflation of the edifice (Figure 8; Supplementary Figure A3), deflation signals are evident, and there is no further increase in volcano-tectonic earthquakes. Vents or dykes are sealed by solidification and cooling, and solidification and crystallization occur within magma reservoirs. When no magma is supplied from the depth to the shallow magma system, the volcano remains dormant, with no surface deformation or seismic signals observed (Figure 8; Supplementary Figure A4). After a period of dormancy, the cycle recommences.

      Schematic illustration of eruptive cycle model of two cases: (A) observable eruption and (B) stealthy eruption. The cycle consists of stages 1) pre-eruptive, when shallow magma inflates as it is fed by the deeper source through the conduit, 2) co-eruptive, when eruption starts as tensile stress exceeds the critical threshold (>0N/m2) and deflation follows as a result of release of magma from chamber, 3) post-eruptive, the cooling stage with deflation of magma chamber, and 4) repose, when flux stopped from deeper magma reservoir and the shallow magma chamber crystallize. In the (A) observable eruption case, upward ground deformation due to the inflation of magma chamber can be detected by InSAR/GPS and increase in number of volcano-tectonic earthquakes caused by the Coulomb failure can be observed (stage 1), and deflation ground deformation signal can be detected since the eruption occurs (stage 2). In contrast, in the (B) stealthy eruption case, no ground deformation and earthquake precursors can be detected at both pre-eruptive stage (stage 1) and co-eruptive stage (stage 2) because the displacement and Coulomb failure are too low. Inspired by the figure in Chaussard et al. (2013).

      Ideally, the detection of both earthquakes and ground deformation enables successful forecasts of volcanic eruptions, including detailed predictions of their timing, location, and magnitude. Notable examples include the 1980 eruption of Mount St. Helens, with accurate forecasts of all subsequent eruptions from April 1981 to December 1982 based on seismic and deformation data, leading to precise predictions without false alarms (Malone et al., 1983; Swanson et al., 1983), the 1991 eruption of Mount Pinatubo, where forecasting significantly mitigated hazards and saved thousands of lives (Voight et al., 1998; 1999), and the accurate short-term prediction of the 2000 Hekla eruption (Soosalu et al., 2005), and these cases highlight the importance of integrating seismic, deformation, and other geophysical data for effective eruption forecasting (Segall, 2013). It is worth noting that an unerupted unrest scenario, marked by observable deformation and seismic signals without a subsequent eruption, can arise, e.g., the Westdahl volcano in the Aleutian Islands displayed significant inflation signals, typically considered precursors to eruptions, yet no eruption occurred within the observed timeframes (Lu et al., 2000; Gong et al., 2015).

      In contrast, the stealthy eruption presents a starkly different scenario compared to its observable counterpart, characterized by minimal detectable signals. During the pre-eruptive stage of a stealthy eruption (depicted in Figure 8; Supplementary Figure B1), the wall rock is stressed due to magma ejection from deep sources. The accumulation of volcano-tectonic (VT) earthquakes may not be readily apparent because, according to our models, there is only minimal shear failure around the magma chamber, with no shear failure detected around the conduit. While edifice inflation occurs, the signals of Line of Sight (LOS) displacement are too subtle to be detected with current monitoring techniques like InSAR. In such cases, the magma reservoir pressurizes, and the tensile stress exceeds the strength of the surrounding crust, potentially leading to an eruption without observable surface inflation or an increase in seismic activity. In the co-eruptive stage, despite the opening of a dyke or conduit resulting in the release of gas and magma from the reservoir, there remains no detectable ground deflation or significant seismic activity, highlighting the eruption’s stealthy nature (Figure 8; Supplementary Figure B2). Similar to observable eruptions, this is followed by the post-eruptive (cooling) stage, where vents collapse or close, and the magma in the upper-crustal reservoir solidifies and crystallizes, cooling the reservoir without observable signals of these processes (Figure 8; Supplementary Figure B3). And finally, the repose stage heralds the onset of a new cycle, with magma replenishment from deeper sources leading to the swelling of magma reservoirs and setting the stage for the next eruption cycle (Figure 8; Supplementary Figure B4).

      The stealthy eruption cycle often evades early detection due to its subtle manifestations, posing challenges in monitoring and forecasting with current technological capabilities. The magma system can either be long-lived or transient. To produce stealthy eruptions, a lack of seismicity requires a large reservoir size, more oblate shape, and low flux, while low displacement requires a small reservoir size, more prolate shape, low flux, and a deeper chamber. Considering the host rock rheology, long-lived, thermally primed magma systems allow for a wider range of the parameter space to result in feasible scenarios for stealthy eruption as opposed to transient systems (Figure 7).

      5.2.2 Veniaminof: open or closed and stealthy?

      Many volcanoes that erupt without observed deformation have long been classified as open-system volcanoes (e.g., Chaussard et al., 2013). In open systems, triggers of eruption include minor magma intrusions, conduit over-pressurization, or lava dome destabilization. One interpretation suggests the deeply rooted plumbing systems pressurize before eruption such as Cleveland and Pavlof volcanoes in Alaska (Lu and Dzurisin, 2014), Galeras in Columbia (Fournier et al., 2010), and Stromboli in Italy, an open-conduit volcano with persistent activity >1,000 years fed by deep seated gas-rich magma (Barberi et al., 2009). Systems with established shallow reservoirs, exhibiting low precursory deformation and the existence of persistent or semi-persistent pathways for magma ascent, are exemplified by volcanoes such as Popocatépetl and Colima in Mexico, and Merapi in Indonesia (Chaussard et al., 2013). Like Veniaminof, the 1999 sub-Plinian basaltic eruption of Shishaldin, Alaska also eluded prediction efforts due to the lack of precursory signals of both InSAR-detectable deformation and VT earthquakes (Lu and Dzurisin, 2014; Moran et al., 2006), although later studies suggested that a swarm of deep long-period (LP) earthquakes, coupled with short-period earthquakes, should have been recognized as precursory signals (Power et al., 2004; Rasmussen et al., 2018). Shishaldin’s lack of precursory signals aligns with both closed-system stealthy eruptions and open-system eruptions; yet vapor saturation pressure studies indicate the inclusion compositions are more consistent with open-system degassing, characterized by a predominant release of CO2 before H2O (Rasmussen et al., 2018).

      Should Veniaminof be considered stealthy, fed by a closed magmatic system? Some geophysical observations (lack of precursory seismicity and observations of ground deformation) and the frequent eruptions at Veniaminof may point to an open volcanic system. According to the volcanic activity summary by AVO (Orr et al., 2024), the majority of historical eruptions since 1830 have likely originated from cone A. This apparent recurrence of eruptions from the same vent could suggest the presence of an open conduit system. Loewen et al. (2021) characterize the 2018 eruption as an exemplar of a frequently active open-vent system, which explains several of Veniaminof’s characteristics: nearly continuous, low-level background seismicity punctuated by two periods of elevated tremors linked to increased tachylite textures that indicate the intense ash production and possibly reflecting a deeper conduit to initiate seismic tremors; explosive activity alongside lava flows integrating into the ice cap; and a transition from Strombolian to Hawaiian eruption styles, indicative of varying gas release and magma dynamics. Nonetheless, while seismic tremors and mineral textures indirectly hint at continuous gas emissions, there is no direct evidence of steady outgassing of volcanic gases (SO2, CO2, H2O) through eruptions or even quiescent periods (Rose et al., 2013; Lyons et al., 2011) and a lack of robust documentation on exceptionally high gas emissions (a hallmark of open-vent systems) or detailed studies correlating continuous gas flux with seismic activity during past eruptions. Additionally, the continuous seismic activity indicating magma movement within the volcanic conduit, typical of open-vent volcanism, is not prominent in Veniaminof’s history, as shown by the low number of detected shallow VT earthquakes. The presence of a persistent lava lake or exposed fresh lava in the form of a dome with frequent gas emissions common in some open-vent volcanoes (Rose et al., 2013) remains unconfirmed. Furthermore, detailed investigation of variations in seismic velocity structure shows decreases in the seismic velocities prior to Veniaminof’s 2004 and 2013 eruptions indicating the potential buildup of magmatic fluids in a closed magmatic system (Bennington et al., 2018). The detailed architecture and dynamics of the magma plumbing system, particularly regarding conduit dynamics essential for sustaining an open-vent paradigm, remain unclear. Detailed analyses of infrasound data to predict open-vent eruptions, used at Villarrica, Chile (Johnson et al., 2018), have yet to be applied to Veniaminof. These gaps complicate classifying Veniaminof’s magmatic system and determining how closely it conforms to the open-vent model.

      The pattern of eruption frequency at Veniaminof does not strictly align with the open system hypothesis. Its historical explosive eruptions including those recent events in 2013, 1983, and 1956 (up to VEI 3, Global Volcanism Program, 2021), may result from new magma injections into a shallow chamber, similar to Askja, Iceland, Krakatau, Indonesia, Nevado del Ruiz, Colombia, and St. Helens, United States which are characterized as closed systems (Colucci and Papale, 2021). Furthermore, concerning its absence of long-term, edifice-wide ground deformation, often linked to open-vent systems, it is critical to note that InSAR data for Veniaminof excludes the summit caldera due to ice coverage but the periphery area. Our models suggest that scenarios producing 15–18 mm of LOS displacement at the caldera’s edge, which falls below Sentinel-1 detection threshold, could indicate 40–270 mm deformation at the center, depending on the chamber’s shape, size, depth, and magma flux. Should caldera data become accessible, we might observe ground deformation indicative of stealthy closed-system eruptions (instead of open-vent eruptions), characterized by negligible peripheral deformation and a marginal increase in seismic activity, driven by low magma influx or warm rheology.

      5.3 Implications for future monitoring of volcanic systems

      Recent technological advancements in data collection and analysis have significantly enhanced our capabilities on volcano monitoring and forecasting. Yet, we are still facing complexities in volcano monitoring and forecasting, including the need for approaches that forecast not only the likelihood of an eruption but also its location, magnitude, style, duration, and potential for ash plumes that impact a long distance (Segall, 2013; Acocella, 2014; Bebbington and Jenkins, 2019). The challenges in short-term volcanic eruption forecasting are significantly heightened by the absence of reliable eruptive precursors despite of the expansion of instrumental networks, the complexity of non-linear pre-eruptive behaviors, and the inherent unpredictability of volcanic systems due to unknown parameters and the potential for sudden changes (Sparks et al., 2012; Acocella, 2014). Specifically, data processing and magma system modeling in forecasting efforts still face significant challenges in quantitatively linking monitoring data to the probabilities of future volcanic events, a limited understanding of volcanic physics, geometries and material properties for physics based modeling, specialization of models for mathematical and computational viability that complicate the task of choosing appropriate models for an as-yet unobserved eruption (Marzocchi and Bebbington, 2012; Poland and Anderson, 2020). Although novel machine learning approaches have facilitated forecasting by detecting critical changes in volcanic activity patterns that indicate transitions at the onset or end of volcanic activity (e.g., classification method by Manley et al., 2021) and despite of advancements in algorithms, the application of machine learning faces challenges including the opaqueness of “black box” models, lack of insights into underlying physical mechanisms, and data limitations due to scarce well-monitored eruptions (Reichstein et al., 2019).

      Observable eruptions exhibiting seismicity, deformation, and gas emissions, as depicted in Figure 8A are the subjects of extensive research. As they allow for direct monitoring of precursory signals, there have been successful forecasts of such volcanic eruptions (e.g., Segall, 2013). In cases where eruptions lack detectable precursory seismic activity, alternative indicators such as deformation and increased gas emissions (such as exemplified by Redoubt Volcano) can provide warnings (Roman and Cashman, 2018), and usually average duration of deformation and degassing (932 days and 282 days, respectively) are longer than seismic unrest (197-day) (Phillipson et al., 2013). Additionally, continuous local stress field monitoring, detectable through changes in seismic velocity and shear-wave splitting, may reveal aseismic staging. Meanwhile, small-scale earthquake swarms, indicative of shallow magma chamber recharging, suggest a years-long eruption risk eruption potential. The deep seismic activity may serve as a near-term harbinger of imminent eruption (Roman and Cashman, 2018).

      In scenarios where ground deformation remains undetected by GNSS and InSAR despite sufficient coverage of the volcanic edifice, it is important to recognize these techniques’ inherent limitations, which have been well documented (Arrowsmith et al., 2021). For example, during the episodes of intense volcanic activity at Mount Etna during spectacular lava fountains, no discernible signals were observed in GNSS or InSAR datasets (Currenti and Bonaccorso, 2019; Carleo et al., 2023; Cardone et al., 2024). These techniques often fail to capture ground deformation induced by small, short-term volcanic events with millimeter-to centimeter-scale accuracy, prompting the deployment of high-precision borehole instruments, such as tiltmeters and strainmeters, on volcano flanks to complement GNSS data and better resolve pre-eruptive deformation (Dzurisin, 2007; Linde and Sacks, 1995; Currenti and Bonaccorso, 2019; Carleo et al., 2023; Cardone et al., 2024). Furthermore, continuing technological research is exploring innovative methods for volcano deformation monitoring, including the fiber optic sensing, which is highly sensitive to external physical processes and observables at specific locations (Jousset et al., 2025).

      In addition to ground-deformation techniques, infrasound monitoring can capture signals of volcanic activity that InSAR or GPS fail to detect, as inaudible low-frequency sound waves are generated by processes such as effusive eruptions and lava lake agitation near volcanic vents (Matoza and Fee, 2018). Real-time analysis of these persistent infrasound signals, emblematic of degassing activities within open-vent systems, alongside seismic and gas compositional data, is crucial for providing key insights into degassing styles and persistent degassing pulses for early eruption forecasting, exemplified by the studies on sustained infrasonic tremor observed at volcanoes like Stromboli, Kilauea, Villarrica, and Shishaldin (Johnson et al., 2018). Open-vent volcanoes are also subject to geochemical monitoring of emissions-whether as plumes, fumaroles, through soil diffusion, or via springs (Edmonds, 2021), with measurements like SO2 emissions being particularly informative (Wilkes et al., 2023; Coppola et al., 2019). Moreover, the characteristic “excess degassing” behavior of open-vent volcanoes, which release more gas than contained in the erupted magma, points to a process of endogenous or cryptic growth magmatic system-the intrusion of unerupted magma within the volcano that does not result in substantial long-term surface deformation. This process is detectable through satellite-derived thermal anomalies in addition to gas emission measurements (Coppola et al., 2019).

      5.3.1 Monitoring enhancements at veniaminof volcano

      Enhancing future monitoring and forecasting capabilities for Veniaminof Volcano necessitates overcoming significant challenges posed by the existing technological limitations. The high-elevation, steep-slope terrain and snow-covered landscape of Veniaminof present formidable obstacles to detecting volcanic deformation via InSAR, primarily due to atmospheric noise interference and loss of coherence over the volcano summit. Possible strategies to enhance detection accuracy include enhancing InSAR coherence using long-wavelength SARs with much shorter temporal repeats and employing variogram modeling for precise simulation of residual atmospheric noise (Beker et al., 2023). Moreover, supplementing InSAR with campaign GPS systems in the vicinity of the caldera could provide detailed temporal changes in the volcanic system from pre-eruptive stages to inform future eruption predictions (previous campaign GPS at Veniaminof, e.g., Fournier and Freymueller, 2008; Drooff and Freymueller, 2021). Furthermore, deploying high-precision borehole instruments, such as tiltmeters and strainmeters, on volcano flanks could improve the detection of subtle ground deformation (e.g., Currenti and Bonaccorso, 2019; Cardone et al., 2024). Additionally, AVO has been upgrading the Aleutian seismic network to broadband seismometers with real-time digital transmission since 2003, with the upgrade at Veniaminof completed in 2022 (Power et al., 2020). By applying advanced machine learning methods to seismic signal processing, the ability to detect volcanic precursor signals amidst background noise may be greatly improved (Zhu and Beroza, 2018; Malfante et al., 2018; Bueno et al., 2020).

      For stealthy volcanoes like Veniaminof and others in Alaska (e.g., Loewen et al., 2021) and around the world, numerical modeling is an invaluable tool. Numerical models of magma reservoir evolution significantly aid in deciphering the processes that lead to the accumulation and transport of magma and eruption-triggering mechanisms, essential for forecasting the future behavior of volcanoes (Caricchi et al., 2021). The 2D FEM models in our study, for example, offer insights into the key factors controlling stress accumulation in the host rock of magma reservoirs including the size, shape depth of the magma chamber, and magma supply rate. Thus, they enhance our understanding of the physical and thermal evolution of volcanic plumbing systems and the conditions that lead to rock rupture and eruption. These are critical for forecasting eruption repose times and guiding strategic data collection efforts.

      6 Conclusion

      Our comprehensive study on Veniaminof Volcano offers pivotal insights into the behavior and characteristics of “stealthy” volcanic eruptions: characterized by the absence of detectable seismic or geodetic precursors, despite the occurrence of an eruption. The numerical modeling results of Veniaminof Volcano, exemplifying stealthy volcanic eruptions within closed magmatic systems, illuminate key factors necessary for such eruptions that occur with minimal ground deformation and seismicity rate changes, based on a series of constraints on various parameter combinations, including magma chamber size, shape, depth, and magma flux rate, under both temperature-dependent and non-temperature-dependent host rock rheology.

      Two critical conditions are a low magma flux rate and a warm rheology of the host rock. A reduced magma ascent rate ensures minimal stress on surrounding rock, limiting ground deformation and volcano-tectonic earthquakes. Warm rheology allows gradual magma movement and deformation within detectable thresholds, generating fewer seismic signals.

      Veniaminof’s eruptive cycle, including pre-eruptive, co-eruptive, post-eruptive, and repose stages, provides a framework applicable to other volcanic systems globally. Insights from Veniaminof’s stealthy eruption guide future forecasts and highlight the need to integrate multidisciplinary data and numerical modeling for accurate predictions. This enhances risk mitigation strategies, reducing volcanic hazards' impact on communities. This study serves as a model for similar volcanoes and advances future monitoring and forecasting efforts in volcanology.

      Data availability statement

      The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author.

      Author contributions

      YL: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Validation, Visualization, Writing – original draft, Writing – review and editing. PG: Conceptualization, Funding acquisition, Methodology, Project administration, Supervision, Writing – original draft, Writing – review and editing. ZL: Conceptualization, Funding acquisition, Writing – original draft, Writing – review and editing. JW: Writing – original draft, Writing – review and editing.

      Funding

      The author(s) declare that financial support was received for the research and/or publication of this article. Investigations of volcanic unrest in the Aleutians is supported by grants from the U.S. National Science Foundation (EAR 1752477, EAR 2122745 – Gregg and Li) and NASA (80-NSSC19K-0357 – Gregg, Wang, and Lu).

      The authors thank the anonymous reviewers and V. Acocella whose comments and discussions greatly improved the manuscript. We also thank R. Maguire, M. Loewan, M. Head, and C. Lundstrom for helpful discussions that greatly enhanced this investigation.

      Conflict of interest

      The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

      Generative AI statement

      The author(s) declare that no Generative AI was used in the creation of this manuscript.

      Publisher’s note

      All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

      Supplementary material

      The Supplementary Material for this article can be found online at: /articles/10.3389/feart.2025.1535083/full#supplementary-material

      References Acocella V. (2014). Great challenges in volcanology: how does the volcano factory work? Front. Earth Sci. 2. 10.3389/feart.2014.00004 Acocella V. Ripepe M. Rivalta E. Peltier A. Galetto F. Joseph E. (2023). Towards scientific forecasting of magmatic eruptions. Nat. Rev. Earth Environ. 5, 522. 10.1038/s43017-023-00492-z Arrowsmith J. R. Brodsky E. E. Cooper C. M. Elliott J. L. Fee D. Fischer K. M. (2021). Recommendations for enabling earth science through NSF’s geophysical facility – a portfolio review of EAR seismology and geodesy instrumentation. Alexandria, VA, United States: US National Science Foundation. Bacon C. R. Sisson T. W. Mazdab F. K. (2007). Young cumulate complex beneath Veniaminof caldera, Aleutian arc, dated by zircon in erupted plutonic blocks. Geol. 35, 491. 10.1130/G23446A.1 Barberi F. Civetta L. Rosi M. Scandone R. (2009). Chronology of the 2007 eruption of Stromboli and the activity of the scientific synthesis group. J. Volcanol. Geotherm. Res. 182, 123130. 10.1016/j.jvolgeores.2008.09.019 Bebbington M. S. Jenkins S. F. (2019). Intra-eruption forecasting. Bull. Volcanol. 81, 34. 10.1007/s00445-019-1294-9 Beker T. Ansari H. Montazeri S. Song Q. Zhu X. X. (2023). Deep learning for subtle volcanic deformation detection with InSAR data in central volcanic zone. IEEE Trans. Geosci. Remote Sens. 61, 120. 10.1109/TGRS.2023.3318469 Bennington N. Haney M. Thurber C. Zeng X. (2018). Inferring magma dynamics at Veniaminof Volcano via application of ambient noise. Geophys. Res. Lett. 45. 10.1029/2018GL079909 Biggs J. Ebmeier S. K. Aspinall W. P. Lu Z. Pritchard M. E. Sparks R. S. J. (2014). Global link between deformation and volcanic eruption quantified by satellite imagery. Nat. Commun. 5, 3471. 10.1038/ncomms4471 Bonafede M. Ferrari C. (2009). Analytical models of deformation and residual gravity changes due to a Mogi source in a viscoelastic medium. Tectonophysics 471, 413. 10.1016/j.tecto.2008.10.006 Bueno A. Zuccarello L. Díaz-Moreno A. Woollam J. Titos M. Benítez C. (2020). PICOSS: Python interface for the classification of seismic signals. Comput. and Geosciences 142, 104531. 10.1016/j.cageo.2020.104531 Cameron C. E. Prejean S. G. Coombs M. L. Wallace K. L. Power J. A. Roman D. C. (2018). Alaska volcano observatory alert and forecasting timeliness: 1989–2017. Front. Earth Sci. 6, 86. 10.3389/feart.2018.00086 Cardone M. Cannata A. Giuffrida M. Iozzia A. Minio V. Viccaro M. (2024). Changing magma recharge/discharge dynamics during the 2020–22 lava fountaining activity at Mt. Etna revealed by tilt deformation and volcanic tremor. J. Volcanol. Geotherm. Res. 449 (April), 108074. 10.1016/j.jvolgeores.2024.108074 Caricchi L. Townsend M. Rivalta E. Namiki A. (2021). The build-up and triggers of volcanic eruptions. Nat. Rev. Earth and Environ. 2, 458476. 10.1038/s43017-021-00174-8 Carleo L. Currenti G. Bonaccorso A. (2023). Clustering of eruptive events from high-precision strain signals recorded during the 2020–2022 lava fountains at the Etna volcano (Italy). Nat. Hazards Earth Syst. Sci. 23 (5), 17431754. 10.5194/nhess-23-1743-2023 Cashman K. V. Sparks R. S. J. Blundy J. D. (2017). Vertically extensive and unstable magmatic systems: a unified view of igneous processes. Science 355, eaag3055. 10.1126/science.aag3055 Chaussard E. Amelung F. Aoki Y. (2013). Characterization of open and closed volcanic systems in Indonesia and Mexico using InSAR time series. J. Geophys. Res. Solid Earth 118, 39573969. 10.1002/jgrb.50288 Colucci S. Papale P. (2021). Deep magma transport control on the size and evolution of explosive volcanic eruptions. Front. Earth Sci. 9, 681083. 10.3389/feart.2021.681083 Cooper G. F. Morgan D. J. Wilson C. J. N. (2017). Rapid assembly and rejuvenation of a large silicic magmatic system: insights from mineral diffusive profiles in the Kidnappers and Rocky Hill deposits, New Zealand. Earth Planet. Sci. Lett. 473, 113. 10.1016/j.epsl.2017.05.036 Coppola D. Laiolo M. Massimetti F. Cigolini C. (2019). Monitoring endogenous growth of open-vent volcanoes by balancing thermal and SO2 emissions data derived from space. Sci. Rep. 9, 9394. 10.1038/s41598-019-45753-4 Currenti G. Bonaccorso A. (2019). Cyclic magma recharge pulses detected by high-precision strainmeter data: the case of 2017 inter-eruptive activity at Etna volcano. Sci. Rep. 9 (1), 75537557. 10.1038/s41598-019-44066-w Del Negro C. Currenti G. Scandura D. (2009). Temperature-dependent viscoelastic modeling of ground deformation: application to Etna volcano during the 1993–1997 inflation period. Phys. Earth Planet. Interiors 172, 299309. 10.1016/j.pepi.2008.10.019 Drooff C. Freymueller J. T. (2021). New constraints on slip deficit on the Aleutian megathrust and inflation at Mt. Veniaminof, Alaska from repeat GPS measurements. Geophys. Res. Lett. 48, e2020GL091787. 10.1029/2020GL091787 Dzurisin D. (2003). A comprehensive approach to monitoring volcano deformation as a window on the eruption cycle. Rev. Geophys. 41, 2001RG000107. 10.1029/2001RG000107 Dzurisin D. (2007). Volcano deformation: geodetic monitoring techniques. Berlin: Springer. 10.1007/978-3-540-49302-0 Edmonds M. (2021). “Geochemical monitoring of volcanoes and the mitigation of volcanic gas hazards,” in Forecasting and planning for volcanic hazards, risks, and disasters (Elsevier), 117151. 10.1016/B978-0-12-818082-2.00004-4 Fournier T. Freymueller J. (2008). Inflation detected at Mount Veniaminof, Alaska, with campaign GPS. Geophys. Res. Lett. 35, 2008GL035503. 10.1029/2008GL035503 Ferretti A. Prati C. Rocca F. (2000). Nonlinear subsidence rate estimation using permanent scatterers in differential SAR interferometry. IEEE. Transact. Geosci. Remo. Sensi. 38 (5), 22022212. Fournier T. J. Pritchard M. E. Riddick S. N. (2010). Duration, magnitude, and frequency of subaerial volcano deformation events: new results from Latin America using InSAR and a global synthesis. Geochem. Geophys. Geosystems 11, 2009GC002558. 10.1029/2009GC002558 Global Volcanism Program (2021). “Report on Veniaminof (United States),” in Weekly volcanic activity report, 7–13 July 2021. Editor Sennert S. (Washington, DC, United States: Smithsonian Institution and US Geological Survey). Gong W. Meyer F. J. Lee C.-W. Lu Z. Freymueller J. (2015). Measurement and interpretation of subtle deformation signals at Unimak Island from 2003 to 2010 using weather model-assisted time series InSAR. J. Geophys. Res. Solid Earth 120, 11751194. 10.1002/2014JB011384 Grapenthin R. Freymueller J. T. Serovetnikov S. S. (2013). Surface deformation of Bezymianny Volcano, Kamchatka, recorded by GPS: the eruptions from 2005 to 2010 and long-term, long-wavelength subsidence. J. Volcanol. Geotherm. Res. 263, 5874. 10.1016/j.jvolgeores.2012.11.012 Gregg P. M. De Silva S. L. Grosfils E. B. Parmigiani J. P. (2012). Catastrophic caldera-forming eruptions: thermomechanics and implications for eruption triggering and maximum caldera dimensions on Earth. J. Volcanol. Geotherm. Res. 241–242, 112. 10.1016/j.jvolgeores.2012.06.009 Gregg P. M. Le Mével H. Zhan Y. Dufek J. Geist D. Chadwick W. W. (2018). Stress triggering of the 2005 eruption of sierra negra volcano, galápagos. Geophys. Res. Lett. 45. 10.1029/2018GL080393 Grosfils E. B. McGovern P. J. Gregg P. M. Galgana G. A. Hurwitz D. M. Long S. M. (2015). Elastic models of magma reservoir mechanics: a key tool for investigating planetary volcanism. Spec. Publ. 401, 239267. 10.1144/SP401.2 Gualda G. A. R. Pamukcu A. S. Ghiorso M. S. Anderson A. T. Sutton S. R. Rivers M. L. (2012). Timescales of quartz crystallization and the longevity of the bishop giant magma body. PLoS ONE 7, e37492. 10.1371/journal.pone.0037492 Hanssen R. F. (2001). “Radar system theory and interferometric processing,” in Radar Interferometry: data interpretation and error analysis (Netherlands: Springer), 960. 10.1007/0-306-47633-9_2 Hickey J. Gottsmann J. Del Potro R. (2013). The large-scale surface uplift in the Altiplano-Puna region of Bolivia: a parametric study of source characteristics and crustal rheology using finite element analysis. Geochem. Geophys. Geosystems 14, 540555. 10.1002/ggge.20057 Hooper A. Zebker H. Segall P. Kampes B. (2004). A new method for measuring deformation on volcanoes and other natural terrains using InSAR persistent scatterers. Geophys. Res. Lett. 31, L23611. 10.1029/2004GL021737 Hooper A. (2006). Persistent scatterer radar interferometry for crustal deformation studies and modeling of volcano deformation. PhD thesis. Stanford (California): Stanford University. Hooper A. Segall P. Zebker H. A. (2007). Persistent scatterer radar interferometry for crustal deformation studies and modeling of volcanic deformation. J. Geophys. Res. Solid. Earth. 112, B07407. 10.1029/2006JB004714 Hooper A. Segall P. Zebker H. A. (2007). Persistent scatterer interferometric synthetic aperture radar for crustal deformation analysis, with application to Volcán Alcedo, Galápagos. J. Geophys. Res. Solid. Earth. 112, B07407. 10.1029/2006JB004763 Hooper A. (2008). A multi-temporal InSAR method incorporating both persistent scatterer and small baseline approaches. Geophys Res. Lett. 35, L16302. 10.1029/2008GL034654 Jellinek A. M. DePaolo D. J. (2003). A model for the origin of large silicic magma chambers: precursors of caldera-forming eruptions. Bull. Volcanol. 65, 363381. 10.1007/s00445-003-0277-y Johnson J. B. Watson L. M. Palma J. L. Dunham E. M. Anderson J. F. (2018). Forecasting the eruption of an open-vent volcano using resonant infrasound tones. Geophys. Res. Lett. 45, 22132220. 10.1002/2017GL076506 Jousset P. Currenti G. Murphy S. Eibl E. P. S. Caudron C. Retailleau L. (2025). “Fiber optic sensing for volcano monitoring and imaging volcanic processes,” in Modern Volcano monitoring (Springer). Karakas O. Degruyter W. Bachmann O. Dufek J. (2017). Lifetime and size of shallow magma bodies controlled by crustal-scale magmatism. Nat. Geosci. 10, 446450. 10.1038/ngeo2959 Le Mével H. Gregg P. M. Feigl K. L. (2016). Magma injection into a long-lived reservoir to explain geodetically measured uplift: Application to the 2007–2014 unrest episode at Laguna del Maule volcanic field, Chile. J. Geophys. Res. Solid Earth 121, 60926108. 10.1002/2016JB013066 Linde A. T. Sacks S. I. (1995). “Continuous monitoring of volcanoes with borehole strainmeters,”. Mauna loa revealed: structure, composition, hystory, and hazards. Editors Rhodes J. M. Lockwood J. P. (Washington DC: American Geophysical Union), 171, 171185. 10.1029/gm092p0171 Loewen M. W. Dietterich H. R. Graham N. Izbekov P. (2021). Evolution in eruptive style of the 2018 eruption of Veniaminof volcano, Alaska, reflected in groundmass textures and remote sensing. Bull. Volcanol. 83, 72. 10.1007/s00445-021-01489-6 Lu Z. Dzurisin D. (2014). InSAR imaging of aleutian volcanoes: monitoring a volcanic arc from space. Springer Berlin Heidelberg. 10.1007/978-3-642-00348-6 Lu Z. Wicks C. Dzurisin D. Thatcher W. Freymueller J. T. McNutt S. R. (2000). Aseismic inflation of Westdahl Volcano, Alaska, revealed by satellite radar interferometry. Geophys. Res. Lett. 27, 15671570. 10.1029/1999GL011283 Lyons J. Rose B. Escobar R. (2011). Volcanic hazards and remote sensing in Pacific Latin America. EoS Trans. 92, 196. 10.1029/2011EO230007 Malfante M. Dalla Mura M. Metaxian J.-P. Mars J. I. Macedo O. Inza A. (2018). Machine learning for volcano-seismic signals: challenges and perspectives. IEEE Signal Process. Mag. 35, 2030. 10.1109/MSP.2017.2779166 Malone S. D. Boyko C. Weaver C. S. (1983). Seismic precursors to the Mount St. Helens eruptions in 1981 and 1982. Science 221, 13761378. 10.1126/science.221.4618.1376 Manley G. F. Mather T. A. Pyle D. M. Clifton D. A. Rodgers M. Thompson G. (2021). Machine learning approaches to identifying changes in eruptive state using multi-parameter datasets from the 2006 eruption of Augustine Volcano, Alaska. J. Geophys. Res. Solid Earth 126, e2021JB022323. 10.1029/2021JB022323 Marzocchi W. Bebbington M. S. (2012). Probabilistic eruption forecasting at short and long time scales. Bull. Volcanol. 74, 17771805. 10.1007/s00445-012-0633-x Matoza R. S. Fee D. (2018). The inaudible rumble of volcanic eruptions. Acoust. Today 14 (1). Miller T. P. McGimsey R. G. Richter D. H. Riehle J. R. Nye C. J. Yount M. E. (1998). Catalog of the historically active volcanoes of Alaska (No. 98–582). Open-File Rep. 10.3133/ofr98582 Miller T. P. Smith R. L. (1987). Late Quaternary caldera-forming eruptions in the eastern Aleutian arc, Alaska. Geology 15, 434438. 10.1130/0091-7613(1987)15<434:LQCEIT>2.0.CO;2 Mittal T. Richards M. A. (2019). Volatile degassing from magma chambers as a control on volcanic eruptions. J. Geophys. Res. Solid Earth 124, 78697901. 10.1029/2018JB016983 Moran S. C. Kwoun O. Masterlark T. Lu Z. (2006). On the absence of InSAR-detected volcano deformation spanning the 1995–1996 and 1999 eruptions of Shishaldin Volcano, Alaska. J. Volcanol. Geotherm. Res. 150, 119131. 10.1016/j.jvolgeores.2005.07.013 Newman A. V. Dixon T. H. Ofoegbu G. I. Dixon J. E. (2001). Geodetic and seismic constraints on recent activity at Long Valley Caldera, California: evidence for viscoelastic rheology. J. Volcanol. Geotherm. Res. 105, 183206. 10.1016/S0377-0273(00)00255-9 Orr T. R. Neal C. A. Power J. A. (2024). Observations of the 2021 eruption of Mount Veniaminof. Alaska: U.S. Geological Survey Scientific Investigations Report, 20245014. 10.3133/sir20245014 Phillipson G. Sobradelo R. Gottsmann J. (2013). Global volcanic unrest in the 21st century: an analysis of the first decade. J. Volcanol. Geotherm. Res. 264, 183196. 10.1016/j.jvolgeores.2013.08.004 Poland M. P. Anderson K. R. (2020). Partly cloudy with a chance of lava flows: forecasting volcanic eruptions in the twenty-first century. J. Geophys. Res. Solid Earth 125, e2018JB016974. 10.1029/2018JB016974 Power J. Friberg P. A. Haney M. M. Parker T. Stihler S. D. Dixon J. P. (2019). A unified catalog of earthquake hypocenters and magnitudes at volcanoes in Alaska—1989 to 2018 (No. 2019–5037). Scientific Investigations Report, U.S. Geological Survey. 10.3133/sir20195037 Power J. A. Haney M. M. Botnick S. M. Dixon J. P. Fee D. Kaufman A. M. (2020). Goals and development of the Alaska Volcano Observatory seismic network and application to forecasting and detecting volcanic eruptions. Seismol. Res. Lett. 91, 647659. 10.1785/0220190216 Power J. A. Stihler S. D. White R. A. Moran S. C. (2004). Observations of deep long-period (DLP) seismic events beneath Aleutian arc volcanoes; 1989–2002. J. Volcanol. Geotherm. Res. 138, 243266. 10.1016/j.jvolgeores.2004.07.005 Rasmussen D. J. Plank T. A. Roman D. C. Power J. A. Bodnar R. J. Hauri E. H. (2018). When does eruption run-up begin? Multidisciplinary insight from the 1999 eruption of Shishaldin volcano. Earth Planet. Sci. Lett. 486, 114. 10.1016/j.epsl.2018.01.001 Reichstein M. Camps-Valls G. Stevens B. Jung M. Denzler J. Carvalhais N. Prabhat F. (2019). Deep learning and process understanding for data-driven Earth system science. Nature 566, 195204. 10.1038/s41586-019-0912-1 Roman D. C. Cashman K. V. (2018). Top-down precursory volcanic seismicity: implications for ‘stealth’ magma ascent and long-term eruption forecasting. Front. Earth Sci. 6, 124. 10.3389/feart.2018.00124 Rose W. I. Palma J. L. Delgado Granados H. Varley N. (2013). “Open-vent volcanism and related hazards: overview,” in Understanding open-vent volcanism and related hazards. Editors Rose W. I. Palma J. L. Delgado Granados H. Varley N. (Boulder, CO, United States: Geological Society of America). 10.1130/2013.249800 498, viiXiii. Segall P. (2013). Volcano deformation and eruption forecasting. Spec. Publ. 380, 85106. 10.1144/SP380.4 Sparks R. S. J. Biggs J. Neuberg J. W. (2012). Monitoring volcanoes. Science 335, 13101311. 10.1126/science.1219485 Soosalu H. Einarsson P. Þorbjarnardóttir B. S. (2005). Seismic activity related to the 2000 eruption of the Hekla volcano, Iceland. Bull. Volcanol. 38, 2136. 10.1007/s00445-005-0417-7 Swanson D. A. Casadevall T. J. Dzurisin D. Malone S. D. Newhall C. G. Weaver C. S. (1983). Predicting eruptions at Mount St. Helens, june 1980 through december 1982. Science 221, 13691376. 10.1126/science.221.4618.1369 Turner S. Costa F. (2007). Measuring timescales of magmatic evolution. Elements 3, 267272. 10.2113/gselements.3.4.267 Voight B. Hoblitt R. P. Clarke A. B. Lockhart A. B. Miller A. D. Lynch L. (1998). Remarkable cyclic ground deformation monitored in real-time on Montserrat, and its use in eruption forecasting. Geophys. Res. Lett. 25, 34053408. 10.1029/98GL01160 Voight B. Sparks R. S. J. Miller A. D. Stewart R. C. Hoblitt R. P. Clarke A. (1999). Magma flow instability and cyclic activity at soufriere hills volcano, Montserrat, British west indies. Science 283, 11381142. 10.1126/science.283.5405.1138 Wallace K. L. Cameron C. E. Waythomas C. F. Perreault J. M. (2020). Major-oxide and trace-element geochemical data from the tephra collected on Veniaminof volcano, Alaska (No. RDF 2020-17). Fairbanks, AK, United States: Alaska Division of Geological and Geophysical Surveys. 10.14509/30578 Waythomas C. F. (2021). Simultaneous effusive and explosive cinder cone eruptions at Veniaminof Volcano, Alaska. Volcanica 4, 295307. 10.30909/vol.04.02.295307 Waythomas C. F. Cameron C. E. Miller T. P. (2022). Eruptive history and hazards of Mount Veniaminof, Alaska. U.S. Geological Survey Scientific Investigations Report, 20225075. 10.3133/sir20225075 Wilkes T. C. Pering T. D. Aguilera F. Layana S. Nadeau P. Kern C. (2023). A new permanent, low-cost, low-power SO2 camera for continuous measurement of volcanic emissions. Front. Earth Sci. 11, 1088992. 10.3389/feart.2023.1088992 Zhan Y. Gregg P. M. (2019). How accurately can we model magma reservoir failure with uncertainties in host rock rheology? J. Geophys. Res. Solid Earth 124, 80308042. 10.1029/2019JB018178 Zhu W. Beroza G. C. (2018). PhaseNet: a deep-neural-network-based seismic arrival time picking method. Geophys. J. Int. 10.1093/gji/ggy423
      ‘Oh, my dear Thomas, you haven’t heard the terrible news then?’ she said. ‘I thought you would be sure to have seen it placarded somewhere. Alice went straight to her room, and I haven’t seen her since, though I repeatedly knocked at the door, which she has locked on the inside, and I’m sure it’s most unnatural of her not to let her own mother comfort her. It all happened in a moment: I have always said those great motor-cars shouldn’t be allowed to career about the streets, especially when they are all paved with cobbles as they are at Easton Haven, which are{331} so slippery when it’s wet. He slipped, and it went over him in a moment.’ My thanks were few and awkward, for there still hung to the missive a basting thread, and it was as warm as a nestling bird. I bent low--everybody was emotional in those days--kissed the fragrant thing, thrust it into my bosom, and blushed worse than Camille. "What, the Corner House victim? Is that really a fact?" "My dear child, I don't look upon it in that light at all. The child gave our picturesque friend a certain distinction--'My husband is dead, and this is my only child,' and all that sort of thing. It pays in society." leave them on the steps of a foundling asylum in order to insure [See larger version] Interoffice guff says you're planning definite moves on your own, J. O., and against some opposition. Is the Colonel so poor or so grasping—or what? Albert could not speak, for he felt as if his brains and teeth were rattling about inside his head. The rest of[Pg 188] the family hunched together by the door, the boys gaping idiotically, the girls in tears. "Now you're married." The host was called in, and unlocked a drawer in which they were deposited. The galleyman, with visible reluctance, arrayed himself in the garments, and he was observed to shudder more than once during the investiture of the dead man's apparel. HoME香京julia种子在线播放 ENTER NUMBET 0016www.jnlpsw.com.cn
      www.eopfwj.com.cn
      www.eekpls.com.cn
      www.v55rz.net.cn
      utogwb.com.cn
      neurub.com.cn
      www.okgnd.com.cn
      www.sdrctz.com.cn
      www.wqliyj.com.cn
      inbp.com.cn
      处女被大鸡巴操 强奸乱伦小说图片 俄罗斯美女爱爱图 调教强奸学生 亚洲女的穴 夜来香图片大全 美女性强奸电影 手机版色中阁 男性人体艺术素描图 16p成人 欧美性爱360 电影区 亚洲电影 欧美电影 经典三级 偷拍自拍 动漫电影 乱伦电影 变态另类 全部电 类似狠狠鲁的网站 黑吊操白逼图片 韩国黄片种子下载 操逼逼逼逼逼 人妻 小说 p 偷拍10幼女自慰 极品淫水很多 黄色做i爱 日本女人人体电影快播看 大福国小 我爱肏屄美女 mmcrwcom 欧美多人性交图片 肥臀乱伦老头舔阴帝 d09a4343000019c5 西欧人体艺术b xxoo激情短片 未成年人的 插泰国人夭图片 第770弾み1 24p 日本美女性 交动态 eee色播 yantasythunder 操无毛少女屄 亚洲图片你懂的女人 鸡巴插姨娘 特级黄 色大片播 左耳影音先锋 冢本友希全集 日本人体艺术绿色 我爱被舔逼 内射 幼 美阴图 喷水妹子高潮迭起 和后妈 操逼 美女吞鸡巴 鸭个自慰 中国女裸名单 操逼肥臀出水换妻 色站裸体义术 中国行上的漏毛美女叫什么 亚洲妹性交图 欧美美女人裸体人艺照 成人色妹妹直播 WWW_JXCT_COM r日本女人性淫乱 大胆人艺体艺图片 女同接吻av 碰碰哥免费自拍打炮 艳舞写真duppid1 88电影街拍视频 日本自拍做爱qvod 实拍美女性爱组图 少女高清av 浙江真实乱伦迅雷 台湾luanlunxiaoshuo 洛克王国宠物排行榜 皇瑟电影yy频道大全 红孩儿连连看 阴毛摄影 大胆美女写真人体艺术摄影 和风骚三个媳妇在家做爱 性爱办公室高清 18p2p木耳 大波撸影音 大鸡巴插嫩穴小说 一剧不超两个黑人 阿姨诱惑我快播 幼香阁千叶县小学生 少女妇女被狗强奸 曰人体妹妹 十二岁性感幼女 超级乱伦qvod 97爱蜜桃ccc336 日本淫妇阴液 av海量资源999 凤凰影视成仁 辰溪四中艳照门照片 先锋模特裸体展示影片 成人片免费看 自拍百度云 肥白老妇女 女爱人体图片 妈妈一女穴 星野美夏 日本少女dachidu 妹子私处人体图片 yinmindahuitang 舔无毛逼影片快播 田莹疑的裸体照片 三级电影影音先锋02222 妻子被外国老头操 观月雏乃泥鳅 韩国成人偷拍自拍图片 强奸5一9岁幼女小说 汤姆影院av图片 妹妹人艺体图 美女大驱 和女友做爱图片自拍p 绫川まどか在线先锋 那么嫩的逼很少见了 小女孩做爱 处女好逼连连看图图 性感美女在家做爱 近距离抽插骚逼逼 黑屌肏金毛屄 日韩av美少女 看喝尿尿小姐日逼色色色网图片 欧美肛交新视频 美女吃逼逼 av30线上免费 伊人在线三级经典 新视觉影院t6090影院 最新淫色电影网址 天龙影院远古手机版 搞老太影院 插进美女的大屁股里 私人影院加盟费用 www258dd 求一部电影里面有一个二猛哥 深肛交 日本萌妹子人体艺术写真图片 插入屄眼 美女的木奶 中文字幕黄色网址影视先锋 九号女神裸 和骚人妻偷情 和潘晓婷做爱 国模大尺度蜜桃 欧美大逼50p 西西人体成人 李宗瑞继母做爱原图物处理 nianhuawang 男鸡巴的视屏 � 97免费色伦电影 好色网成人 大姨子先锋 淫荡巨乳美女教师妈妈 性nuexiaoshuo WWW36YYYCOM 长春继续给力进屋就操小女儿套干破内射对白淫荡 农夫激情社区 日韩无码bt 欧美美女手掰嫩穴图片 日本援交偷拍自拍 入侵者日本在线播放 亚洲白虎偷拍自拍 常州高见泽日屄 寂寞少妇自卫视频 人体露逼图片 多毛外国老太 变态乱轮手机在线 淫荡妈妈和儿子操逼 伦理片大奶少女 看片神器最新登入地址sqvheqi345com账号群 麻美学姐无头 圣诞老人射小妞和强奸小妞动话片 亚洲AV女老师 先锋影音欧美成人资源 33344iucoom zV天堂电影网 宾馆美女打炮视频 色五月丁香五月magnet 嫂子淫乱小说 张歆艺的老公 吃奶男人视频在线播放 欧美色图男女乱伦 avtt2014ccvom 性插色欲香影院 青青草撸死你青青草 99热久久第一时间 激情套图卡通动漫 幼女裸聊做爱口交 日本女人被强奸乱伦 草榴社区快播 2kkk正在播放兽骑 啊不要人家小穴都湿了 www猎奇影视 A片www245vvcomwwwchnrwhmhzcn 搜索宜春院av wwwsee78co 逼奶鸡巴插 好吊日AV在线视频19gancom 熟女伦乱图片小说 日本免费av无码片在线开苞 鲁大妈撸到爆 裸聊官网 德国熟女xxx 新不夜城论坛首页手机 女虐男网址 男女做爱视频华为网盘 激情午夜天亚洲色图 内裤哥mangent 吉沢明歩制服丝袜WWWHHH710COM 屌逼在线试看 人体艺体阿娇艳照 推荐一个可以免费看片的网站如果被QQ拦截请复制链接在其它浏览器打开xxxyyy5comintr2a2cb551573a2b2e 欧美360精品粉红鲍鱼 教师调教第一页 聚美屋精品图 中韩淫乱群交 俄罗斯撸撸片 把鸡巴插进小姨子的阴道 干干AV成人网 aolasoohpnbcn www84ytom 高清大量潮喷www27dyycom 宝贝开心成人 freefronvideos人母 嫩穴成人网gggg29com 逼着舅妈给我口交肛交彩漫画 欧美色色aV88wwwgangguanscom 老太太操逼自拍视频 777亚洲手机在线播放 有没有夫妻3p小说 色列漫画淫女 午间色站导航 欧美成人处女色大图 童颜巨乳亚洲综合 桃色性欲草 色眯眯射逼 无码中文字幕塞外青楼这是一个 狂日美女老师人妻 爱碰网官网 亚洲图片雅蠛蝶 快播35怎么搜片 2000XXXX电影 新谷露性家庭影院 深深候dvd播放 幼齿用英语怎么说 不雅伦理无需播放器 国外淫荡图片 国外网站幼幼嫩网址 成年人就去色色视频快播 我鲁日日鲁老老老我爱 caoshaonvbi 人体艺术avav 性感性色导航 韩国黄色哥来嫖网站 成人网站美逼 淫荡熟妇自拍 欧美色惰图片 北京空姐透明照 狼堡免费av视频 www776eom 亚洲无码av欧美天堂网男人天堂 欧美激情爆操 a片kk266co 色尼姑成人极速在线视频 国语家庭系列 蒋雯雯 越南伦理 色CC伦理影院手机版 99jbbcom 大鸡巴舅妈 国产偷拍自拍淫荡对话视频 少妇春梦射精 开心激动网 自拍偷牌成人 色桃隐 撸狗网性交视频 淫荡的三位老师 伦理电影wwwqiuxia6commqiuxia6com 怡春院分站 丝袜超短裙露脸迅雷下载 色制服电影院 97超碰好吊色男人 yy6080理论在线宅男日韩福利大全 大嫂丝袜 500人群交手机在线 5sav 偷拍熟女吧 口述我和妹妹的欲望 50p电脑版 wwwavtttcon 3p3com 伦理无码片在线看 欧美成人电影图片岛国性爱伦理电影 先锋影音AV成人欧美 我爱好色 淫电影网 WWW19MMCOM 玛丽罗斯3d同人动画h在线看 动漫女孩裸体 超级丝袜美腿乱伦 1919gogo欣赏 大色逼淫色 www就是撸 激情文学网好骚 A级黄片免费 xedd5com 国内的b是黑的 快播美国成年人片黄 av高跟丝袜视频 上原保奈美巨乳女教师在线观看 校园春色都市激情fefegancom 偷窥自拍XXOO 搜索看马操美女 人本女优视频 日日吧淫淫 人妻巨乳影院 美国女子性爱学校 大肥屁股重口味 啪啪啪啊啊啊不要 操碰 japanfreevideoshome国产 亚州淫荡老熟女人体 伦奸毛片免费在线看 天天影视se 樱桃做爱视频 亚卅av在线视频 x奸小说下载 亚洲色图图片在线 217av天堂网 东方在线撸撸-百度 幼幼丝袜集 灰姑娘的姐姐 青青草在线视频观看对华 86papa路con 亚洲1AV 综合图片2区亚洲 美国美女大逼电影 010插插av成人网站 www色comwww821kxwcom 播乐子成人网免费视频在线观看 大炮撸在线影院 ,www4KkKcom 野花鲁最近30部 wwwCC213wapwww2233ww2download 三客优最新地址 母亲让儿子爽的无码视频 全国黄色片子 欧美色图美国十次 超碰在线直播 性感妖娆操 亚洲肉感熟女色图 a片A毛片管看视频 8vaa褋芯屑 333kk 川岛和津实视频 在线母子乱伦对白 妹妹肥逼五月 亚洲美女自拍 老婆在我面前小说 韩国空姐堪比情趣内衣 干小姐综合 淫妻色五月 添骚穴 WM62COM 23456影视播放器 成人午夜剧场 尼姑福利网 AV区亚洲AV欧美AV512qucomwwwc5508com 经典欧美骚妇 震动棒露出 日韩丝袜美臀巨乳在线 av无限吧看 就去干少妇 色艺无间正面是哪集 校园春色我和老师做爱 漫画夜色 天海丽白色吊带 黄色淫荡性虐小说 午夜高清播放器 文20岁女性荫道口图片 热国产热无码热有码 2015小明发布看看算你色 百度云播影视 美女肏屄屄乱轮小说 家族舔阴AV影片 邪恶在线av有码 父女之交 关于处女破处的三级片 极品护士91在线 欧美虐待女人视频的网站 享受老太太的丝袜 aaazhibuo 8dfvodcom成人 真实自拍足交 群交男女猛插逼 妓女爱爱动态 lin35com是什么网站 abp159 亚洲色图偷拍自拍乱伦熟女抠逼自慰 朝国三级篇 淫三国幻想 免费的av小电影网站 日本阿v视频免费按摩师 av750c0m 黄色片操一下 巨乳少女车震在线观看 操逼 免费 囗述情感一乱伦岳母和女婿 WWW_FAMITSU_COM 偷拍中国少妇在公车被操视频 花也真衣论理电影 大鸡鸡插p洞 新片欧美十八岁美少 进击的巨人神thunderftp 西方美女15p 深圳哪里易找到老女人玩视频 在线成人有声小说 365rrr 女尿图片 我和淫荡的小姨做爱 � 做爱技术体照 淫妇性爱 大学生私拍b 第四射狠狠射小说 色中色成人av社区 和小姨子乱伦肛交 wwwppp62com 俄罗斯巨乳人体艺术 骚逼阿娇 汤芳人体图片大胆 大胆人体艺术bb私处 性感大胸骚货 哪个网站幼女的片多 日本美女本子把 色 五月天 婷婷 快播 美女 美穴艺术 色百合电影导航 大鸡巴用力 孙悟空操美少女战士 狠狠撸美女手掰穴图片 古代女子与兽类交 沙耶香套图 激情成人网区 暴风影音av播放 动漫女孩怎么插第3个 mmmpp44 黑木麻衣无码ed2k 淫荡学姐少妇 乱伦操少女屄 高中性爱故事 骚妹妹爱爱图网 韩国模特剪长发 大鸡巴把我逼日了 中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片中国张柏芝做爱片 大胆女人下体艺术图片 789sss 影音先锋在线国内情侣野外性事自拍普通话对白 群撸图库 闪现君打阿乐 ady 小说 插入表妹嫩穴小说 推荐成人资源 网络播放器 成人台 149大胆人体艺术 大屌图片 骚美女成人av 春暖花开春色性吧 女亭婷五月 我上了同桌的姐姐 恋夜秀场主播自慰视频 yzppp 屄茎 操屄女图 美女鲍鱼大特写 淫乱的日本人妻山口玲子 偷拍射精图 性感美女人体艺木图片 种马小说完本 免费电影院 骑士福利导航导航网站 骚老婆足交 国产性爱一级电影 欧美免费成人花花性都 欧美大肥妞性爱视频 家庭乱伦网站快播 偷拍自拍国产毛片 金发美女也用大吊来开包 缔D杏那 yentiyishu人体艺术ytys WWWUUKKMCOM 女人露奶 � 苍井空露逼 老荡妇高跟丝袜足交 偷偷和女友的朋友做爱迅雷 做爱七十二尺 朱丹人体合成 麻腾由纪妃 帅哥撸播种子图 鸡巴插逼动态图片 羙国十次啦中文 WWW137AVCOM 神斗片欧美版华语 有气质女人人休艺术 由美老师放屁电影 欧美女人肉肏图片 白虎种子快播 国产自拍90后女孩 美女在床上疯狂嫩b 饭岛爱最后之作 幼幼强奸摸奶 色97成人动漫 两性性爱打鸡巴插逼 新视觉影院4080青苹果影院 嗯好爽插死我了 阴口艺术照 李宗瑞电影qvod38 爆操舅母 亚洲色图七七影院 被大鸡巴操菊花 怡红院肿么了 成人极品影院删除 欧美性爱大图色图强奸乱 欧美女子与狗随便性交 苍井空的bt种子无码 熟女乱伦长篇小说 大色虫 兽交幼女影音先锋播放 44aad be0ca93900121f9b 先锋天耗ばさ无码 欧毛毛女三级黄色片图 干女人黑木耳照 日本美女少妇嫩逼人体艺术 sesechangchang 色屄屄网 久久撸app下载 色图色噜 美女鸡巴大奶 好吊日在线视频在线观看 透明丝袜脚偷拍自拍 中山怡红院菜单 wcwwwcom下载 骑嫂子 亚洲大色妣 成人故事365ahnet 丝袜家庭教mp4 幼交肛交 妹妹撸撸大妈 日本毛爽 caoprom超碰在email 关于中国古代偷窥的黄片 第一会所老熟女下载 wwwhuangsecome 狼人干综合新地址HD播放 变态儿子强奸乱伦图 强奸电影名字 2wwwer37com 日本毛片基地一亚洲AVmzddcxcn 暗黑圣经仙桃影院 37tpcocn 持月真由xfplay 好吊日在线视频三级网 我爱背入李丽珍 电影师傅床戏在线观看 96插妹妹sexsex88com 豪放家庭在线播放 桃花宝典极夜著豆瓜网 安卓系统播放神器 美美网丝袜诱惑 人人干全免费视频xulawyercn av无插件一本道 全国色五月 操逼电影小说网 good在线wwwyuyuelvcom www18avmmd 撸波波影视无插件 伊人幼女成人电影 会看射的图片 小明插看看 全裸美女扒开粉嫩b 国人自拍性交网站 萝莉白丝足交本子 七草ちとせ巨乳视频 摇摇晃晃的成人电影 兰桂坊成社人区小说www68kqcom 舔阴论坛 久撸客一撸客色国内外成人激情在线 明星门 欧美大胆嫩肉穴爽大片 www牛逼插 性吧星云 少妇性奴的屁眼 人体艺术大胆mscbaidu1imgcn 最新久久色色成人版 l女同在线 小泽玛利亚高潮图片搜索 女性裸b图 肛交bt种子 最热门有声小说 人间添春色 春色猜谜字 樱井莉亚钢管舞视频 小泽玛利亚直美6p 能用的h网 还能看的h网 bl动漫h网 开心五月激 东京热401 男色女色第四色酒色网 怎么下载黄色小说 黄色小说小栽 和谐图城 乐乐影院 色哥导航 特色导航 依依社区 爱窝窝在线 色狼谷成人 91porn 包要你射电影 色色3A丝袜 丝袜妹妹淫网 爱色导航(荐) 好男人激情影院 坏哥哥 第七色 色久久 人格分裂 急先锋 撸撸射中文网 第一会所综合社区 91影院老师机 东方成人激情 怼莪影院吹潮 老鸭窝伊人无码不卡无码一本道 av女柳晶电影 91天生爱风流作品 深爱激情小说私房婷婷网 擼奶av 567pao 里番3d一家人野外 上原在线电影 水岛津实透明丝袜 1314酒色 网旧网俺也去 0855影院 在线无码私人影院 搜索 国产自拍 神马dy888午夜伦理达达兔 农民工黄晓婷 日韩裸体黑丝御姐 屈臣氏的燕窝面膜怎么样つぼみ晶エリーの早漏チ○ポ强化合宿 老熟女人性视频 影音先锋 三上悠亚ol 妹妹影院福利片 hhhhhhhhsxo 午夜天堂热的国产 强奸剧场 全裸香蕉视频无码 亚欧伦理视频 秋霞为什么给封了 日本在线视频空天使 日韩成人aⅴ在线 日本日屌日屄导航视频 在线福利视频 日本推油无码av magnet 在线免费视频 樱井梨吮东 日本一本道在线无码DVD 日本性感诱惑美女做爱阴道流水视频 日本一级av 汤姆avtom在线视频 台湾佬中文娱乐线20 阿v播播下载 橙色影院 奴隶少女护士cg视频 汤姆在线影院无码 偷拍宾馆 业面紧急生级访问 色和尚有线 厕所偷拍一族 av女l 公交色狼优酷视频 裸体视频AV 人与兽肉肉网 董美香ol 花井美纱链接 magnet 西瓜影音 亚洲 自拍 日韩女优欧美激情偷拍自拍 亚洲成年人免费视频 荷兰免费成人电影 深喉呕吐XXⅩX 操石榴在线视频 天天色成人免费视频 314hu四虎 涩久免费视频在线观看 成人电影迅雷下载 能看见整个奶子的香蕉影院 水菜丽百度影音 gwaz079百度云 噜死你们资源站 主播走光视频合集迅雷下载 thumbzilla jappen 精品Av 古川伊织star598在线 假面女皇vip在线视频播放 国产自拍迷情校园 啪啪啪公寓漫画 日本阿AV 黄色手机电影 欧美在线Av影院 华裔电击女神91在线 亚洲欧美专区 1日本1000部免费视频 开放90后 波多野结衣 东方 影院av 页面升级紧急访问每天正常更新 4438Xchengeren 老炮色 a k福利电影 色欲影视色天天视频 高老庄aV 259LUXU-683 magnet 手机在线电影 国产区 欧美激情人人操网 国产 偷拍 直播 日韩 国内外激情在线视频网给 站长统计一本道人妻 光棍影院被封 紫竹铃取汁 ftp 狂插空姐嫩 xfplay 丈夫面前 穿靴子伪街 XXOO视频在线免费 大香蕉道久在线播放 电棒漏电嗨过头 充气娃能看下毛和洞吗 夫妻牲交 福利云点墦 yukun瑟妃 疯狂交换女友 国产自拍26页 腐女资源 百度云 日本DVD高清无码视频 偷拍,自拍AV伦理电影 A片小视频福利站。 大奶肥婆自拍偷拍图片 交配伊甸园 超碰在线视频自拍偷拍国产 小热巴91大神 rctd 045 类似于A片 超美大奶大学生美女直播被男友操 男友问 你的衣服怎么脱掉的 亚洲女与黑人群交视频一 在线黄涩 木内美保步兵番号 鸡巴插入欧美美女的b舒服 激情在线国产自拍日韩欧美 国语福利小视频在线观看 作爱小视颍 潮喷合集丝袜无码mp4 做爱的无码高清视频 牛牛精品 伊aⅤ在线观看 savk12 哥哥搞在线播放 在线电一本道影 一级谍片 250pp亚洲情艺中心,88 欧美一本道九色在线一 wwwseavbacom色av吧 cos美女在线 欧美17,18ⅹⅹⅹ视频 自拍嫩逼 小电影在线观看网站 筱田优 贼 水电工 5358x视频 日本69式视频有码 b雪福利导航 韩国女主播19tvclub在线 操逼清晰视频 丝袜美女国产视频网址导航 水菜丽颜射房间 台湾妹中文娱乐网 风吟岛视频 口交 伦理 日本熟妇色五十路免费视频 A级片互舔 川村真矢Av在线观看 亚洲日韩av 色和尚国产自拍 sea8 mp4 aV天堂2018手机在线 免费版国产偷拍a在线播放 狠狠 婷婷 丁香 小视频福利在线观看平台 思妍白衣小仙女被邻居强上 萝莉自拍有水 4484新视觉 永久发布页 977成人影视在线观看 小清新影院在线观 小鸟酱后丝后入百度云 旋风魅影四级 香蕉影院小黄片免费看 性爱直播磁力链接 小骚逼第一色影院 性交流的视频 小雪小视频bd 小视频TV禁看视频 迷奸AV在线看 nba直播 任你在干线 汤姆影院在线视频国产 624u在线播放 成人 一级a做爰片就在线看狐狸视频 小香蕉AV视频 www182、com 腿模简小育 学生做爱视频 秘密搜查官 快播 成人福利网午夜 一级黄色夫妻录像片 直接看的gav久久播放器 国产自拍400首页 sm老爹影院 谁知道隔壁老王网址在线 综合网 123西瓜影音 米奇丁香 人人澡人人漠大学生 色久悠 夜色视频你今天寂寞了吗? 菲菲影视城美国 被抄的影院 变态另类 欧美 成人 国产偷拍自拍在线小说 不用下载安装就能看的吃男人鸡巴视频 插屄视频 大贯杏里播放 wwwhhh50 233若菜奈央 伦理片天海翼秘密搜查官 大香蕉在线万色屋视频 那种漫画小说你懂的 祥仔电影合集一区 那里可以看澳门皇冠酒店a片 色自啪 亚洲aV电影天堂 谷露影院ar toupaizaixian sexbj。com 毕业生 zaixian mianfei 朝桐光视频 成人短视频在线直接观看 陈美霖 沈阳音乐学院 导航女 www26yjjcom 1大尺度视频 开平虐女视频 菅野雪松协和影视在线视频 华人play在线视频bbb 鸡吧操屄视频 多啪啪免费视频 悠草影院 金兰策划网 (969) 橘佑金短视频 国内一极刺激自拍片 日本制服番号大全magnet 成人动漫母系 电脑怎么清理内存 黄色福利1000 dy88午夜 偷拍中学生洗澡磁力链接 花椒相机福利美女视频 站长推荐磁力下载 mp4 三洞轮流插视频 玉兔miki热舞视频 夜生活小视频 爆乳人妖小视频 国内网红主播自拍福利迅雷下载 不用app的裸裸体美女操逼视频 变态SM影片在线观看 草溜影院元气吧 - 百度 - 百度 波推全套视频 国产双飞集合ftp 日本在线AV网 笔国毛片 神马影院女主播是我的邻居 影音资源 激情乱伦电影 799pao 亚洲第一色第一影院 av视频大香蕉 老梁故事汇希斯莱杰 水中人体磁力链接 下载 大香蕉黄片免费看 济南谭崔 避开屏蔽的岛a片 草破福利 要看大鸡巴操小骚逼的人的视频 黑丝少妇影音先锋 欧美巨乳熟女磁力链接 美国黄网站色大全 伦蕉在线久播 极品女厕沟 激情五月bd韩国电影 混血美女自摸和男友激情啪啪自拍诱人呻吟福利视频 人人摸人人妻做人人看 44kknn 娸娸原网 伊人欧美 恋夜影院视频列表安卓青青 57k影院 如果电话亭 avi 插爆骚女精品自拍 青青草在线免费视频1769TV 令人惹火的邻家美眉 影音先锋 真人妹子被捅动态图 男人女人做完爱视频15 表姐合租两人共处一室晚上她竟爬上了我的床 性爱教学视频 北条麻妃bd在线播放版 国产老师和师生 magnet wwwcctv1024 女神自慰 ftp 女同性恋做激情视频 欧美大胆露阴视频 欧美无码影视 好女色在线观看 后入肥臀18p 百度影视屏福利 厕所超碰视频 强奸mp magnet 欧美妹aⅴ免费线上看 2016年妞干网视频 5手机在线福利 超在线最视频 800av:cOm magnet 欧美性爱免播放器在线播放 91大款肥汤的性感美乳90后邻家美眉趴着窗台后入啪啪 秋霞日本毛片网站 cheng ren 在线视频 上原亚衣肛门无码解禁影音先锋 美脚家庭教师在线播放 尤酷伦理片 熟女性生活视频在线观看 欧美av在线播放喷潮 194avav 凤凰AV成人 - 百度 kbb9999 AV片AV在线AV无码 爱爱视频高清免费观看 黄色男女操b视频 观看 18AV清纯视频在线播放平台 成人性爱视频久久操 女性真人生殖系统双性人视频 下身插入b射精视频 明星潜规测视频 mp4 免賛a片直播绪 国内 自己 偷拍 在线 国内真实偷拍 手机在线 国产主播户外勾在线 三桥杏奈高清无码迅雷下载 2五福电影院凸凹频频 男主拿鱼打女主,高宝宝 色哥午夜影院 川村まや痴汉 草溜影院费全过程免费 淫小弟影院在线视频 laohantuiche 啪啪啪喷潮XXOO视频 青娱乐成人国产 蓝沢润 一本道 亚洲青涩中文欧美 神马影院线理论 米娅卡莉法的av 在线福利65535 欧美粉色在线 欧美性受群交视频1在线播放 极品喷奶熟妇在线播放 变态另类无码福利影院92 天津小姐被偷拍 磁力下载 台湾三级电髟全部 丝袜美腿偷拍自拍 偷拍女生性行为图 妻子的乱伦 白虎少妇 肏婶骚屄 外国大妈会阴照片 美少女操屄图片 妹妹自慰11p 操老熟女的b 361美女人体 360电影院樱桃 爱色妹妹亚洲色图 性交卖淫姿势高清图片一级 欧美一黑对二白 大色网无毛一线天 射小妹网站 寂寞穴 西西人体模特苍井空 操的大白逼吧 骚穴让我操 拉好友干女朋友3p