TIME-DEPENDENT MODELING OF VOLCANO DEFORMATION IN ALASKA AND TRANSIENT DETECTION USING MACHINE LEARNING METHODS By Xueming Xue A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of Environmental Geosciences—Doctor of Philosophy 2021 ABSTRACT TIME-DEPENDENT MODELING OF VOLCANO DEFORMATION IN ALASKA AND TRANSIENT DETECTION USING MACHINE LEARNING METHODS By Xueming Xue Geodetic observations provide a measure of surface deformation caused by underground geophysical processes. GPS and InSAR are the most two widely used geodetic techniques, and are often combined because of their complementary advantages. This dissertation presents several methods and study cases to detect and model transient deformation recorded by GPS and InSAR measurements. A time-dependent inversion filter based on the unscented Kalman filter is developed to combine InSAR and GPS measurements for volcano deformation modeling. The InSAR and GPS data at Okmok since the 2008 eruption shows a shallow magma system to be sill-like structure at ~0.9 km and a Mogi source at ~3.2 km depth. The shallow sill structure was not seen before the 2008 eruption. Five inflation episodes have been observed from 2008 to 2019, each decaying exponentially in time. A hydraulic model with the episodically increasing pressure of the deep magma reservoir and increasing influx from deep to shallow sources volcano is discussed to explain consecutive inflation events. Most of the available InSAR and GPS data, including campaign GPS, continuous GPS, single InSAR and InSAR times series, are used to recover the magma supply history of Akutan, Makushin and Westdahl volcano over the past 20-25 years since the first geodetic measurements were made there. The time-dependent volume changes show diverse temporal patterns across the three volcanoes. A small residual velocity signal is extracted by comparing the estimated regional velocity to the current block models. Two possible coupling models are tested to explain the residual velocity signal as a result of slip deficit on the megathrust. A sequence-to-sequence ML method is proposed to detect the transient signals in GPS time series. The model is trained using synthetic data time series that are relatively short in length but can be applied to time series of arbitrary length. The model is validated by being applied to some real-world datasets containing transient signals with various characteristics. By decimating a time series to a lower data rate to effectively compress time, the model is able to recognize long-duration transients even when it has been trained only on short-duration transients. This indicates the similarity in nature between the short and long transients. ACKNOWLEDGMENTS Finally, I get to write this part of the dissertation. From 2015 to 2021, I started my journey in Alaska, then moved to Michigan, and witnessed a global pandemic. I am so happy at the moment not only because of this dissertation but also some many things that happened and people that I knew during the past years. Those experiences totally changed my life as well as my view of the world. I am thankful to be a student of Jeff Freymueller. Jeff is a great advisor who is always patient and provides me many opportunities to explore what I am interested in. I also have a great appreciation for my committee of Shawn Wei, Kevin Mackey, Allen McNamara, and Min Chen who all contributed to my study. Thanks also goes to Mike West, Carl Tape, and Jessica Larsen for their guidance back in Alaska. I would like to acknowledge Zhong Lu, Feifei Qu, and Wenyu Gong who provide InSAR data and answer my questions. I would like to thank my fellow students from our geodesy lab, Shanshan Li, Hugh Harper, Connor Drooff. You are awesome people to work together. I learned a lot from you. I am lucky to know and be surrounded by many friends who love so much. There are so many meaningful memories with Benjamin Smith, Sheng Zhao and Kay Wang. Special thanks go to Carl and Cathy Cady for they make me feel I have a family in this country. I wish I could but there is not enough space to list all the names here. Last but most important, I would like to thank my loving wife Mengyan who is the biggest support and encouragement throughout the time. Without you, I would have never made it. You earned this degree along with me. Thanks also goes out to my family in China. You are so far but you always make me feel I am with you. iv TABLE OF CONTENTS LIST OF TABLES ........................................................................................................................ vii LIST OF FIGURES ..................................................................................................................... viii Chapter 1: Introduction ................................................................................................................... 1 Chapter 2: Modeling the post-eruptive deformation at Okmok based on the GPS and InSAR timeseries: changes in the shallow magma storage system ............................................................ 6 2.1 Abstract ................................................................................................................................. 6 2.2 Introduction ........................................................................................................................... 6 2.2.1 Okmok volcano and its recent eruptions ........................................................................ 6 2.2.2 Time-dependent inversion filters ................................................................................... 8 2.3 Data ..................................................................................................................................... 10 2.3.1 GPS data ....................................................................................................................... 10 2.3.2 InSAR data ................................................................................................................... 12 2.4 Time-dependent inversion method ..................................................................................... 14 2.4.1 Volcanic source models ............................................................................................... 14 2.4.2 Nuisance parameters associated with the measurements ............................................. 15 2.4.3 Initialization and hyperparameter tuning ..................................................................... 17 2.4.4 Relative weights between datasets ............................................................................... 18 2.5 Time-dependent source models .......................................................................................... 19 2.5.1 GPS only ...................................................................................................................... 19 2.5.2 InSAR only .................................................................................................................. 21 2.5.3 GPS and InSAR ........................................................................................................... 23 2.6 Volume change from 2008 to 2019 .................................................................................... 27 2.7 Discussion ........................................................................................................................... 28 2.7.1 Subsidence within the caldera region........................................................................... 28 2.7.2 Shallow magma sources ............................................................................................... 31 2.7.3 Magma supply dynamics ............................................................................................. 34 2.7.4 How close is Okmok to the next eruption? .................................................................. 36 2.8 Conclusion .......................................................................................................................... 37 Chapter 3: 25-year history of volcano magma supply in the east-central Aleutian arc, Alaska ... 39 3.1 Abstract ............................................................................................................................... 39 3.2 Introduction ......................................................................................................................... 39 3.3 Data ..................................................................................................................................... 42 3.4 Method ................................................................................................................................ 44 3.5 Results ................................................................................................................................. 46 3.6 Discussion ........................................................................................................................... 51 3.6.1 Magma Supply ............................................................................................................. 51 3.6.2 Plate Motion and Slip Deficit on the Megathrust ........................................................ 54 3.7 Conclusion .......................................................................................................................... 57 v Chapter 4: Single-station detection of transient deformation in GPS time series based on machine learning methods ........................................................................................................................... 60 4.1 Abstract ............................................................................................................................... 60 4.2 Introduction ......................................................................................................................... 60 4.3 Methods .............................................................................................................................. 62 4.3.1 Synthetic data generation and labeling ........................................................................ 62 4.3.2 Model structure, training and prediction ...................................................................... 64 4.4 Results ................................................................................................................................. 66 4.4.1 Testing on synthetic data.............................................................................................. 66 4.4.2 Application to a real dataset ......................................................................................... 68 4.4.3 Testing on additional data that are “unseen” in the training data ................................ 71 4.5 Conclusion .......................................................................................................................... 73 Chapter 5: Conclusion................................................................................................................... 75 BIBLIOGRAPHY ......................................................................................................................... 78 vi LIST OF TABLES Table 2.1. Estimated parameters to construct the covariance matrix for InSAR images. ............ 13 Table 3.1. Estimated regional velocities relative to North America in comparison with block model estimates. ....................................................................................................................................... 51 Table 3.2. GPS site locations. ....................................................................................................... 58 vii LIST OF FIGURES Figure 1.1. Example of GPS and InSAR time series. (left) GPS data at site OKSO located at south flank of Okmok volcano. (right) InSAR data is from TerraSAR descending track 93 for Okmok volcano. ........................................................................................................................................... 1 Figure 2.1. Location and shaded relief map of Okmok volcano. Cone A, Cone Ahmanilix and the GPS sites used in our study are labelled. ........................................................................................ 7 Figure 2.2. (a) Baseline measurements from the three continuous sites OKCE, OKNC and OKSO relative to the site OKFG. (b) Examples of the InSAR data from each track showing line-of-sight displacements. Black arrow on the top-left panel is satellite azimuth and white arrow is satellite look angle. (c) Data availability in time. Each bar represents a data point. The gray line separates the GPS data into two parts. The data from 2008 to 2016 are used in the combination with InSAR mainly to determine the source geometry, and the data after are used only for the volume change estimation (see section 2.4). .......................................................................................................... 11 Figure 2.3. Likelihood (-2L) contours as function of the hyperparameters for volume change and source location (a) when using GPS data only and (b) when using InSAR data only. Star markers represent the preferred model by maximum likelihood estimation. ............................................. 20 Figure 2.4. Time-dependent source model based on the GPS data only. All parameters are allowed to vary with time except for the ratio. The solid lines are the estimated values and the error bars show 2σ uncertainties. The reference point of the Cartesian coordinate is the geometric center of the caldera. The volume change is assumed to be zero in the beginning. .................................... 21 Figure 2.5. Time-dependent source model based on the InSAR data only. The axis ratio is fixed to be 1. The gray lines are the GPS estimates shown in Figure 2.4. The uncertainties of the GPS estimates are small and not shown here. ....................................................................................... 22 Figure 2.6. Time-dependent source model based on the GPS and InSAR data. (a) Estimates of the deep Mogi source. All the parameters are allowed to vary with time except for the depth. (b) Estimates of the shallow sill. All the parameters are estimated as constant except for the volume change. .......................................................................................................................................... 25 Figure 2.7. Deformation caused by the Mogi source at (a) GPS sites, (c) TS93 and (e) T116. Deformation caused by the sill source at (b) GPS sites, (d) TS93 and (f) T116. Horizontal and vertical GPS displacements are shown by black and red vectors, respectively. InSAR displacements are shown in line-of-sight direction. ..................................................................... 26 Figure 2.8. Total volume change from the sill and the Mogi source. The exponential fit and 95% confidence are denoted by solid and dashed lines, respectively. The time scale τ and characteristic strength ΔV of each inflation event is shown in the figure. .......................................................... 28 Figure 2.9. (a) Residuals of the GPS time series. (b) Residuals of selected interferograms in Figure 2. The linear phase ramp of each interferogram has been removed. The estimated initial position and the common mode has been removed for each time series. ................................................... 29 viii Figure 2.10. (a) OKCE monument viewed from Cone E looking NE toward the gates of the caldera and Cone D. The 1997 lava flow is visible on the left. (b) OKNC monument viewed looking NNE toward to Cone B which is visible to the right. The monument is installed in rock, which is covered by ash from the 2008 deposit. (c-f) GPS timeseries projected into LOS direction in comparison with the InSAR LOS timeseries. Random walks at the GPS sites and linear phase ramps have been estimated and removed. ................................................................................................................ 30 Figure 2.11. Geographic and depth map of the sources. The red star represents the Mogi source and the red dashed circle (line) represents the sill source obtained by combining the GPS and InSAR data. The blue circle is the spheroid source given by the GPS data. For comparison, the black circle denotes the pre-eruptive source. The Mogi source determined by the InSAR data only (section 2.5.2) is close to the Mogi source in the joint inversion (red star). ................................. 32 Figure 2.12. Conceptual model of the magma plumbing system at Okmok and dynamics of deep magma sources. The shallow sill source is inferred to have provided some magma (likely the more evolved component) for the 2008 eruption (Freymueller and Kaufman, 2010; Larsen et al., 2013) but had not been connected to the main shallow magma body prior to the 2008 eruption. Based on a hydraulic model (Lengliné et al., 2008), our results indicate that the pressure of the deep source may be episodic and the magma influx from the deep to shallow source may increase over time. ....................................................................................................................................................... 36 Figure 3.1. Tectonic setting of Aleutian arc, Alaska (inset) and topographic map of the study region (black rectangle in inset map) which includes Westdahl, Akutan and Makushin volcanoes. The red lines show the rigid blocks for Alaska Peninsula and eastern Aleutians from Elliott and Freymueller (2020). The brown polygon outlines the 1957 earthquake aftershock zone (House et al., 1981). The focal mechanism represents the 2003 M6.6 Unimak Island earthquake, which happened 200 km east of Dutch harbor just outside the study region. The triangles show the locations of the continuous (red) and campaign (blue) GPS sites. The red vectors denote the average GPS velocity at the continuous sites, with the estimated regional velocity removed. .... 40 Figure 3.2. Time-dependent volume change estimates with 2σ uncertainties (error bars) for (a) Westdahl, (b) Akutan and (c) Makushin. The gray curves or lines are fit to the estimates. Other47 Figure 3.3. Volume change estimates for Makushin volcano in the bestfit model and the models when the east component of the regional velocity is fixed to the values deviated 1 or 3 times of the estimation uncertainty from the original estimate in the paper. The two black boxes denote the time windows that are affected by the parameter tradeoff between the regional velocity and volume change because of a lack of InSAR or near-field GPS data. ........................................................ 49 Figure 3.4. Two alternate plate coupling models: (a) slip deficit only at shallow depth (<15 km) or (b) slip deficit only at 35-50 km. The dashed line outlines the Alaska trench. Slab contours are shown in black lines, and the outline of the locked zone as a red rectangle with the slip deficit rate given as a percentage of plate motion rate. The top left inset denotes the secular velocity from our estimates (red), previous block models (green) and their difference (blue). The percentage number is the slip deficit rate relative to the convergent rate between the Pacific plate and Peninsula block (6.6 cm/yr). The predicted deformation caused by the slip deficit is denoted as black vectors. .. 56 ix Figure 4.1. An example of (a) raw synthetic timeseries, (b) detrended synthetic timeseries and (c) the labels for the timeseries........................................................................................................... 63 Figure 4.2. A diagram showing the structure and workflow of the ML model. Each input dn is a data point of a timeseries. The bidirectional LSTM layer and the dense layer are used to extract features. The final output gives the transient probabilities. .......................................................... 65 Figure 4.3. (top) Precision and (bottom) recall curve as a function of transient size. .................. 67 Figure 4.4. (top) The east component of displacement and predicted transient probabilities for station ALBH on southern Vancouver Island. (bottom) Transient detections for all stations. Shown are the transient probability greater than 0.5 in any of the three components (east, north, up). ... 69 Figure 4.5. Transient velocity field, remodeled based on transient detections with probability greater than 0.5, in the Cascadia from 2005 to 2016 from the results of (left) our ML model, (middle) the RSI model and (right) their difference. .................................................................... 71 Figure 4.6. Test results for AB48, AC03 and ATW2. For AC03, the gray line in (b) is the prediction results using the original time series as input and the yellow line shows the results using the resampled time series. ................................................................................................................... 74 x Chapter 1: Introduction The Global Position System (GPS) and the Interferometric Synthetic Aperture Radar (InSAR) are the two satellite-based techniques for measuring surface displacements that are widely used for geodetic studies. GPS provides three-dimensional (east, north and vertical in a local coordinate system) displacement measurements of in-situ ground monuments with high temporal resolution (usually daily solutions). InSAR can provide displacement measurements with high spatial resolution in the satellite line-of-sight (LOS) direction but at low temporal resolution. An example of GPS time series and InSAR time series are shown in Figure 1.1. Figure 1.1. Example of GPS and InSAR time series. (left) GPS data at site OKSO located at south flank of Okmok volcano. (right) InSAR data is from TerraSAR descending track 93 for Okmok volcano. InSAR time series are used more and more often together with GPS time series for investigating how geophysical processes evolves in time (e.g. Bekaert et al., 2016; Henderson & 1 Pritchard, 2013). However, both GPS and InSAR time series can contain significant nuisance signals that are not of interest (e.g. phase ramps, seasonal variations, common mode errors). There has been a sustaining interest in the geodetic community in how to detect and model transient deformation of interest in GPS and InSAR timeseries measurements, and there are two typical strategies for this question. One way is to directly model the timeseries time-dependently, such as the network inversion filter (NIF) (Segall & Matthews, 1997). The NIF uses a stochastic description to describe how a geophysical process evolves in time using the linear Kalman filter. The extended NIF (ENIF) (McGuire & Segall, 2003), which is essentially based on the extended Kalman filter, was later proposed for nonlinear estimation of the hyperparameters involved in the model, such as smoothing wieghts. The NIF method has achieved great success in studying slow slip events for the Cascadia region (e.g. Bartlow et al., 2014; Schmidt & Gao, 2010). The NIF is a sequential estimator in what the modeling process ingests the observations epoch by epoch, providing an evolving estimate of the geophysical parameters. It is more computational efficient compared to the one-step least squares. In terms of volcano deformation modeling, the unscented Kalman filter (UKF) (Julier & Uhlmann, 2004; Wan & Van Der Merwe, 2000) and the ensemble Kalman filter (EnKF) (Evensen, 1994) have been adopted for nonlinear models (Fournier et al., 2009; Gregg & Pettijohn, 2016). Instead of calculating the exact Jacobians of a nonlinear model, the UKF/EnKF approximates the Jacobians using simulation results after passing a set of points selected randomly or according to a certain rule. Therefore, they perform better than the NIF/ENIF for modeling volcano deformation, as the volcanic models are often highly nonlinear. One focus of this dissertation is to address two main questions using geodetic observations on volcano deformation: what is the geometry of the magma plumbing system beneath a volcano 2 and how does the magma intrusion evolve over time? A different geometry of the magma plumbing system will result in a different deformation pattern in space. The geometry of the magma body or bodies that are being pressurized can be determined with InSAR data and GPS data from multiple sites distributed across the volcano. The temporal evolution of the magma intrusion can be measured by the relative volume change in time, which can be estimated using time series observations. The time-dependent volume change can be then used to study the magma supply dynamics based on more physics-based models. In Chapter 2, I aim to expand the time-dependent inversion filter based on the UKF that was first developed by (Fournier et al., 2009) for modeling volcano deformation. The original method is able to utilize GPS observations only, here it is improved to be able to integrate InSAR observations as well. I then apply the improved inversion filter to modeling the post-eruptive deformation at Okmok volcano, Alaska, since 2008. By using this improved inversion filter, I address discrepancies in the source depth inferred by GPS and InSAR individually (Freymueller & Kaufman, 2010; Qu et al., 2015), and invert for time-dependent volume change estimates from 2008 to 2019. Based on volume change history, the volume deficit theory for forecasting the eruption timing at Okmok and magma influx dynamics from deep to shallow sources are discussed. This work has been published in Journal of Geophysical Research: Solid Earth (doi:10.1029/2019JB017801). Chapter 3 presents another work using the method proposed in Chapter 2, where I utilize most of the available InSAR and GPS data to recover the magma supply history of the three active volcanoes (i.e. Westdahl, Akutan and Makushin) in the esat-centeral Aleutian arc, Alaska, over the 20–25 years since the first geodetic measurements were made at each volcano. This study aims to provide an overall picture of volume change over time by combining GPS and InSAR data that 3 were from different satellite missions or viewing geometries, or focused only on a specific time window. Another goal of this work is to evaluate the subduction coupling effect in the east-central Aleutians. The block motion of the A laska Peninsula and eastern Aleutians is well known (Elliott & Freymueller, 2020; Li et al., 2016), but not to the central part due to limited GPS observations which are also affected by the large amount of transient volcanic deformation. With the volcano deformation is well modeled, an accurate secular velocity for the region is determined and can be used to evaluate the subduction coupling model. This work has been published in Geophysical Research Letters (doi:10.1029/2020GL088388). The alternative way to model time-dependent deformation is to detect and extract the transient deformation first, and then the extracted transient deformation is then used for modeling the physical model of the process involved. Proposed methods include principal component analysis (Dong et al., 2006), Gaussian wavelet transform (Melbourne et al., 2005), sparse estimation with template base functions (Riel et al., 2014), multichannel singular spectrum analysis (Damian Walwer et al., 2016) and relative strength index (RSI) (Crowell et al., 2016). While the former way, based on some kind of Kalman filter method, uses the time series data to estimate parameters of a physical model, these methods are based on a pure mathematical decomposition for transient signal detection, aiming to improve the signal-to-noise ratio by utilizing the spatial or temporal correlation inherent to the time series. More recently, machine learning (ML) methods, especially neural networks, are used as a data mining tool in many geophysical studies, such as seismic phase picking (Ross et al., 2018; Zhu & Beroza, 2018), volcanic deformation detecting in satellite images (Anantrasirichai et al., 2019), and slow slip events detecting in sea floor pressure data (He et al., 2020), etc. 4 Inspired by the above ML application, Chapter 4 presents a sequence-to- sequence ML model for transient detection based on long short-term memory (LSTM) cells, a special type of recurrent neural network (RNN). In the previous ML applications for time series observations (e.g. He et al., 2020; Zhu & Beroza, 2018), the common solution is to scan the time series through a time window and then do the classification just using the data within that time window. This solution seems not an efficient way as the neighboring time windows are using a lot of the same data. Also it cannot take into account global information from the entire time series, which might be more useful for the detection. In this work, I propose a sequence to sequence ML model which scans through the entire time series only once, taking each data point as the input to output its detection result. The methods for model training and parameter optimization are presented. The model performance is evaluated through synthetic tests, and then validated on some real data. Chapter 5 summarizes the major findings of this dissertation and discuss some works could be explored for improving on this work. 5 Chapter 2: Modeling the post-eruptive deformation at Okmok based on the GPS and InSAR timeseries: changes in the shallow magma storage system 2.1. Abstract Based on the unscented Kalman filter (UKF), we develop a time-dependent inversion filter combining GPS and InSAR time series observations for modeling volcano deformation. We use the Variance Component Estimation method as a means to assign the relative weights for GPS and InSAR data. Then we use the inversion filter to model the post-eruptive deformation at Okmok volcano, Alaska. We find that a Mogi source at 3-4 km depth fits the InSAR data well, while the best fit to the GPS data is an oblate spheroid source at about 2.5 km depth. Our final model consists of a shallow sill at ~0.9 km and a Mogi source at ~3.2 km depth, which well fit both the GPS and InSAR data simultaneously. We think the Mogi source obtained here is the same source account for the pre-eruptive deformation. The shallow sill is a new structure that was not seen before the 2008 eruption. From 2008 to 2019, we have observed five inflation episodes, each of which decays exponential in time. We find that the characteristic time scale of those inflation episodes decreases with respect to time. The total volume change from the two sources is 0.068 km3, which recovers 50-60% of the volume decrease during the 2008 eruption. 2.2. Introduction 2.2.1. Okmok volcano and its recent eruptions Okmok volcano is one of the most active volcanoes in the central Aleutian arc, Alaska. Within the past century, significant eruptions occurred in 1945, 1958 and 1997 (Miller et al., 1998). These events were all effusive eruptions originating from Cone A (Figure 2.1). The most recent eruption occurred in 2008 from July 12 to late August. There are many differences between the 2008 eruption and the previous three events. First, the 2008 eruption opened a new cone, 6 Ahmanilix, close to Cone D but far from Cone A. Second, it was more phreatomagmatic rather than effusive, indicating the 2008 eruption was involved with water to a great extent. Several hundred million cubic meters of ash and tephra deposits were ejected during the eruption and blanketed the whole caldera and part of the northeastern flank. Third, the size of the 2008 eruption was bigger than the previous historical eruptions. It produced a minimum dense rock equivalent volume of 0.17 km3 which is about three times larger than the volumes of the 1958 and 1997 eruptions (J. F. Larsen et al., 2015). The geodetic models also show the total magma volume decrease during the 2008 eruption was 2~3 times larger than the volume decrease during the 1997 eruption (Lu et al., 2010). Figure 2.1. Location and shaded relief map of Okmok volcano. Cone A, Cone Ahmanilix and the GPS sites used in our study are labelled. 7 Persistent ground deformation has also been observed between the eruptions at Okmok volcano. (Fournier et al., 2009) inverted 6 years of continuous and campaign GPS data and found a stable Mogi source at ~2.6 km depth. (Lu et al., 2010) analyzed 150 InSAR images between the 1997 and the 2008 eruptions and found a Mogi source at ~3 km depth. (Biggs et al., 2010)analyzed a combined dataset including more than 15 years of InSAR and continuous GPS measurement and found that all pre-2008 data could be explained well by a single Mogi source at 3.4 km depth. These studies, no matter what data are used, indicate that all the deformation before the 2008 eruption can be explained by a single Mogi source centered beneath the caldera at ~3 km depth. As for the post-eruptive source model, (Freymueller & Kaufman, 2010) found that a shallower source than the pre-eruptive source located at 1.9 km depth was best to interpret the co-eruptive and early post-eruptive deformation based on GPS measurements. However, (Qu et al., 2015) inverted the post-eruptive InSAR time series observations from 2008 to 2014 and found a source at ~3.9 km depth, apparently twice as deep as the GPS estimated source. Although these two studies focused on different time windows, the source locations they determined seem to be too different, which prompts a further study to determine the source geometry for the post-eruptive period. 2.2.2. Time-dependent inversion filters It is always better to take advantage of time-series observations, which capture the temporal evolution of the source model and then can help us understand the magma dynamics. The most common way to utilize time series observations is through a time-dependent inversion filter. Compared to least squares methods, inversion filters invert the time series of observations in a sequential way and demand much less computational cost. There have been a lot of efforts made to develop such techniques. The pioneering application in geodesy is the network inversion filter 8 (NIF), which is essentially a linear Kalman filter. The NIF was first proposed to study slow slip sequences and has been successfully applied in many cases (Bartlow et al., 2014; Bekaert et al., 2016; Segall & Matthews, 1997). The extended NIF is then developed to simultaneously estimate the slip distribution and Laplacian smoothing factor based on the extended Kalman filter which linearizes the system by calculating the Jacobians (McGuire & Segall, 2003). Although the extended NIF can deal with nonlinear systems, there will be usually unpredictable errors when linearizing highly nonlinear systems which are exactly the case for most volcanic source models. Apparently, the NIFs are not suitable to invert for volcanic sources because the forward models of volcano deformation are usually highly nonlinear. Another extension of the Kalman filter, called the unscented Kalman filter (UKF), was later introduced to model volcano deformation by (Fournier et al., 2009). They developed an inversion filter that was able to model GPS time series and then recovered the kinemics of the volcanic source at Okmok volcano from 2004 to 2008. Instead of calculating the exact Jacobians, the UKF approximates the Jacobians of the nonlinear model by passing through a set of sigma point selected based on a method called the unscented transformation (Julier & Uhlmann, 2004; Wan & Van Der Merwe, 2000). In this study, we will extend the work of (Fournier et al., 2009) by making the UKF be able to combine GPS and InSAR measurements. Readers should refer to (Fournier et al., 2009) or (Julier & Uhlmann, 2004) for technical details about UKF. (Gregg & Pettijohn, 2016) recently introduced a data assimilation framework for inverting volcanic sources using Finite Element Models (FEMs) and the ensemble Kalman filter (EnKF). Their method has been also demonstrated to be effective studying real datasets (Albright et al., 2019; Zhan & Gregg, 2017). The EnKF does not calculate the Jacobians explicitly either. It approximates the Jacobians of the nonlinear system using an ensemble of models, which are 9 sampled by the use of the Monte Carlo Markov Chain method. Both the UKF and EnKF will perform well when the models are highly nonlinear, but we chose to use UKF in this work because (Fournier et al., 2009) have demonstrated that it is effective for this kind of problem. In addition, the UKF is more efficient when using simple models although both the UKF and EnKF can be expected to perform well and give equivalent results. The number of sampling points used in the UKF is 2n+1, where n is the number of unknown parameters, and thus for our problem with very few parameters the number of samples required is small. When n becomes very large, the UKF becomes less efficient than the EnKF because the number of models needed in the EnKF ensemble is less dependent on the number of unknown parameters. (Gregg & Pettijohn, 2016) recommended using between 50 to 200 models in the ensemble with the EnKF. Also, it requires a significant convergence time although this should decrease as the number of models in the ensemble decreases. We expect that the EnKF should perform better after reaching convergence for FEMs with a larger number of parameters. 2.3. Data 2.3.1. GPS data GPS measurements record the displacement of an in-situ monument in three dimensions (east, north and up). The post-eruptive GPS network at Okmok volcano is composed of seven sites (Figure 2.1). Among them, OKCE, OKSO and OKFG are continuous sites installed in 2002 and OKNC was built in 2010 after the eruption (Figure 2.2a). OK13, OK31 and OK37 are temporary sites surveyed for several months right after the eruption, and were resurveyed in the summer of 2010. Data availability in time is shown in Figure 2.2c. All the GPS data were processed based on daily solutions and transformed into ITRF2008, using the approach of (Fu & Freymueller, 2012). 10 Following (Fournier et al., 2009), we use weekly averaged positions and rescale the noise covariance based on the short-term scatter. The scaling factor is 4.5. Figure 2.2. (a) Baseline measurements from the three continuous sites OKCE, OKNC and OKSO relative to the site OKFG. (b) Examples of the InSAR data from each track showing line-of-sight displacements. Black arrow on the top-left panel is satellite azimuth and white arrow is satellite look angle. (c) Data availability in time. Each bar represents a data point. The gray line separates the GPS data into two parts. The data from 2008 to 2016 are used in the combination with InSAR mainly to determine the source geometry, and the data after are used only for the volume change estimation (see section 2.4). 11 We use the baselines relative to OKFG instead of the absolute displacements. First, the GPS network is not robust enough for the inversion filter to directly estimate regional velocity. OKFG is the only site located far from the caldera center and most of the signals recorded should be non-volcanic deformation. By using the baseline observations, we avoid estimating the regional velocity and can have a better constraint on the volcanic sources. Additionally, deeper deflation were observed at OKFG during the first year after the eruption (Freymueller & Kaufman, 2010). These deflation signals may reflect activities of deeper magma sources. It is not likely that we can have a reliable estimation for that with our data. Using the baseline measurements also remove the signals from potential much deep sources so that we can have stronger constraints on the shallow sources which we will focus on in this paper. Note that the GPS measurement equation (which will be later discussed in section 3.2) when using baseline measurements should be written as dbaseline = 𝑔(𝐦) − 𝑔OKFG (𝐦), where 𝑔 is the forward model of a volcanic source model. The GPS measurements mentioned in the rest of this paper refers to the baseline measurement relative to OKFG. 2.3.2. InSAR data InSAR measurements record the displacement of a coherent pixel on the ground in the line- of-sight direction to a remote sensing satellite. We use the InSAR data from Envisat ascending track 222 (E222), Envisat descending track 115 (E115), TerraSAR ascending track 116 (T116), and TerraSAR descending track 93 (T93) (Figure 2.2b). The Envisat data are available until 2010 and the TerraSAR data are available until 2014 (Figure 2.2c). The data were processed(Qu et al., 2015), and readers are referred to that paper for details on satellite information and data processing. We use their data, but mask out the InSAR coherent pixels inside the caldera (Figure 2.2b) because we find that the InSAR data would be fit more poorly if those pixels are included, and there are 12 systematic post-fit residuals inside the caldera. We infer that these residuals represent actual surficial deformation related to the thick tephra deposits from the 2008 eruption (see section 2.7.1). Table 2.1. Estimated parameters to construct the covariance matrix for InSAR images. Track Name Covariance Length Scale Acquisition Time Variance (mm2) (Reference Time) (mm2) (km) 20090810 23.212 1.095 0.679 20090821 48.431 1.518 0.456 20090901 45.170 2.041 0.559 T93 20110806 151.107 4.104 0.714 (20080914) 20110817 189.402 4.838 0.964 20110908 108.336 2.607 0.679 20140901 354.503 24.856 0.919 20140912 340.063 21.716 0.865 20120917 27.800 2.173 0.687 20121009 34.559 2.338 0.633 20130711 31.879 1.287 0.694 T116 20130904 47.824 2.690 1.115 (20110702) 20140811 47.685 2.111 0.813 20140822 41.869 1.986 1.167 20140924 78.075 4.160 1.224 20141005 42.212 1.740 0.685 20090901 10.566 0.852 0.879 EI115 20100713 7.929 0.730 2.698 (20090623) 20100817 24.868 1.311 2.367 20100921 47.020 6.300 2.746 20090701 37.839 2.515 1.389 EI222 20090805 33.327 2.171 1.589 (20081029) 20091014 31.539 0.810 1.022 13 Unlike GPS, the uncertainties of InSAR data are not directly given and it is not trivial to include InSAR uncertainties in the time-dependent inversion filter. The InSAR uncertainties describe the observation noise which are mainly caused by atmospheric delays and are spatially correlated. They are assumed to follow a Gaussian distribution, with correlations between nearby pixels so that the covariance matrix is fully populated instead of diagonal. Here we use the method of sample semi-variogram and semi-covariogram to construct the variance-covariance matrix for each interferogram (Lohman & Simons, 2005). First, we estimate a Mogi source and a linear phase ramp for each original interferogram using least squares. The residuals should be mainly the errors. We then carry out the error analysis on the residuals. The parameters for constructing the variance- covariance matrix of each interferogram are listed in Table 2.1. Finally, we resample all the interferograms using the quadtree method (Simons et al., 2002) and propagate the corresponding variance-covariance matrix. These downsampled interferograms, with their variance-covariance matrix, are used in any time-dependent inversions of this study. 2.4. Time-dependent inversion method 2.4.1. Volcanic source models In this study, we use the analytical model which calculates the deformation from an arbitrary oriented spheroid (P. F. Cervelli, 2013) but we simplify the model to make the inversions more robust. First, we fix the dip to 90° to only estimate a vertical oriented spheroid source. In this case, the strike becomes a dummy variable. Then, we estimate the ratio of the vertical to the horizontal axis instead of the absolute semi-axes dimensions (a point source assumption), which allows us to know the shape but not the actual size of the source. The point source assumption is valid when the size of the source is small enough compared to its depth. For very shallow sources, 14 the point source assumption would be incorrect. We have tested if the point source assumption is valid when the source is estimated to be shallow enough. The simplified model then determines a vertically oriented spheroid using five parameters: location (east, north, vertical), volume change and axis ratio, This model is equivalent to a Mogi source when the axis ratio equals to 1 (Mogi, 1958). When the ratio is larger than 1, the model is a prolate spheroidal source (Yang et al., 1988). When the ratio is smaller than 1, the model is equivalent to an oblate spheroidal source when the ratio is less than 1. When the ratio is much less than 1 (e.g. 0.01), the model approximates a penny-shaped crack source or so called circular sill (Fialko et al., 2001). We present the source strength in terms of volume change. This is an interpretation that the pressure changes are due to increasing the volume of the reservoir by adding new magma. 2.4.2. Nuisance parameters associated with the measurements While always the same source parameters, different nuisance parameters will be estimated for different types of data. For GPS data, an initial position for each component (east, north and vertical) is estimated along with the volcanic source parameters. As we are using baseline observations, we do not need to estimate a regional velocity or frame translation but we estimate a common mode error for all GPS data. The source of common mode could be the partitioning between OKFG and other sites due to the slab subduction, considering that OKFG is more than 10 km apart and is closer to the Aleutian trench. For InSAR data, two types of nuisance parameters are estimated. One is the linear phase ramp. As the phase ramps should not be correlated in time, they are estimated as white noise. The other type of nuisance parameter is the initial volume change, which represents the volume change at the reference acquisition time of each track. They are estimated as a constant. If more than one 15 source is modeled in the UKF, an initial volume change will be estimated for each source. For the sake of reducing the number of unknown parameters, the initial volume change parameters can be shared by multiple InSAR timeseries if their starting times are same or very close (e.g. shorter than one week). This is not vital but can help reduce the total number of unknown parameters, especially when modeling multiple sources. In this study, we only estimate one set of initial volume change parameter for both T93 and E222 as their reference acquisition times are in the same week. Given one GPS timeseries and one InSAR timeseries, now we can write the state transition equation with the process noise as below, 𝐦k 𝐈 𝐦k−1 Δ𝑡𝑘 𝜹2𝐦 𝐜k 𝟎 𝐜k−1 𝜹2𝐜 𝐝0 = 𝐈 𝐝0 , 𝐐𝐤 = 𝟎 (1) 𝐫k 0 𝐫k−1 𝜹2𝐫 [ 𝑉0 ] [ 1 ] [ 𝑉0 ] [ 0] where 𝐦 is general source model parameters; 𝐜 is the common mode for all GPS timeseries; 𝐝0 is the initial position for each site; 𝐫 is the linear phase ramp and 𝑉0 is the reference volume change for the InSAR timeseries; 𝐈 is the identity matrix. The process noise deviation is denoted as 𝜹. The corresponding measurement equation at a certain epoch can be written as below, 𝐝GPS G(𝐦, 𝐱 GPS ) + 𝐜 + 𝐝0 𝜎 𝐆𝐏𝐒 𝝐2GPS [ ]=[ ],𝐑 = [ ] (2) 𝐝InSAR G(𝐦, 𝑉0 , 𝐱 InSAR ) + 𝐋𝐫 𝜎 𝐈𝐧𝐒𝐀𝐑 𝝐2InSAR where 𝐝 is the data; 𝝐 is the data noise deviation; x is the measurement location; 𝐋𝐫 represents a linear phase ramp where 𝐋 = [𝐈 𝐱 𝐈𝐧𝐒𝐀𝐑 𝟏 𝐱 𝐈𝐧𝐒𝐀𝐑 𝟐 ]; G is the forward model. The scaling factor of the noise covariance matrix is 𝜎 𝐆𝐏𝐒 for GPS data and 𝜎 𝐈𝐧𝐒𝐀𝐑 for InSAR data. Other symbols are same as in equation (1). As the time interval of InSAR time series can be large, (Bekaert et al., 2016) introduced pseudo observations when integrating InSAR data into the NIF. In practice, it is also helpful to use 16 pseudo observations for InSAR data in the UKF. The pseudo observations essentially provide additional constraints on the initial/reference volume change and then help calibrate all the InSAR data to the same time frame along with GPS data. (Bekaert et al., 2016) generated the pseudo observations to be all zeros. We find that this may be robust for the linear inversion filters (i.e. the network inversion filter), but may cause numerical instabilities for a nonlinear inversion filter (i.e. the UKF). Here we generate InSAR pseudo observations from a normal distribution with zero mean and very small deviation (e.g. 1e-6). 2.4.3. Initialization and hyperparameter tuning An inappropriate prior or initialization could cause the divergence of the UKF. However, based on extensive testing for this problem, it is not necessary to give a highly accurate prior because a reasonable prior is easy to obtain, and a prior with uncertainty of a few kilometers for the horizontal location and 1-2 km for the source depth ensured rapid convergence. In choosing the sigma points for the UKF, we used hyperparameters determined by (Fournier et al., 2009), which also worked well for the similar source models used in this study, and convergence of the UKF was not sensitive to variations in these hyperparameters. It is much more important to determine the optimal hyperparameters (i.e. 𝜹𝐦 , 𝜹𝐜 and 𝜹𝐫 ) that control the temporal smoothness of the corresponding unknown parameters. A larger hyperparameter allows the parameter to change faster in time but could cause overfitting of noise in the data. A smaller hyperparameter makes the parameter change slower but could cause underfitting of the data. We find that the values of 𝜹𝐜 and 𝜹𝐫 for the nuisance parameters do not impact much the final results, so we simply set 𝜹𝐜 = [10,10,10] (mm) and 𝜹𝐜 = [100,1,1] (mm, mm/km, mm/km). The values of 𝜹𝐦 for the source parameter plays a much more significant role on the final estimates. There are two typical methods used to obtain the optimal hyperparameters: 17 direct estimation (McGuire & Segall, 2003) and maximum likelihood estimation (Segall & Matthews, 1997). Direct estimation is more compact and efficient but requires a denser and more stable observation network to get a reliable result. When combining the GPS and InSAR measurements, the direct estimation method is not very robust because the observation network is dramatically changing whenever an InSAR image is available, unless the time interval of a InSAR timeseries is comparable with GPS (e.g. weekly solution). Therefore, we choose the maximum likelihood estimation to estimate the optimal values of 𝜹𝐦 in our study. 2.4.4. Relative weights between datasets Finally, when combining multiple datasets we apply a variance scaling on each dataset, using the simplified Helmert Variance Component Estimation (HVCE) approach (Xu et al., 2009)[Xu et al., 2009]. From equation (2), we expect the data noise to be statistically described by the noise covariance matrix. In other words, the reduced chi-square of each dataset properly weighted is expected to be equal to 1. The final scaling factor of each dataset are computed using an iterative process, starting with equal scaling factors for all datasets. At each iteration, the scale factor of the data covariance matrix for each dataset is the square root of the reduced chi-square, based on the residuals of the data type. 𝜎𝑗 = 1, 𝑗=1 T 𝜎𝑗 = v𝑗−1 Qv𝑗−1 (3) √ , 𝑗>1 { 𝑛 where 𝑗 counts the iteration number; v is the residual for a given data type; n is the number of data; and Q is the inverse of the noise covariance matrix for that data. The residuals of the 𝑗th iteration, v𝑗−1 were calculated with 𝜎𝑗−1 being the scale factor and then a new scale factor 𝜎𝑗 is calculated based on v𝑗−1. The iteration stops when 𝜎𝑗 = 𝜎𝑗−1. 18 2.5. Time-dependent source models In this section, we present the time-dependent source models that best fit the data in each case (i.e. GPS only, InSAR only, GPS and InSAR jointly). We also describe how to determine the source geometry and the hyperparameters for each case. For the comparison with InSAR, we use the GPS data before 2016 because we do not have InSAR data after then. 2.5.1. GPS only For inverting the GPS data, we use the simplified spheroid model and estimate the axis ratio. Following the general source model for pre-eruptive deformation, we tried to invert a single Mogi source (the axis ratio is fixed to be 1). We found that a single Mogi source is not adequate to fit the data because of two main reasons. (1) With different initializations, either OKCE or OKNC will be underestimated. There are more than 10% of the signals in the vertical component that are not captured, which is a large amount of signal considering the total magnitude of deformation at these two sites. (2) There were unreasonable shifts in the estimates of source location estimates when OKCE temporarily dropped out of the network in 2015, indicating that the OKCE and OKNC data are best fit by two different Mogi sources. The next step is to determine the hyperparameters. We first estimate the axis ratio as a constant, which assumes the source shape does not change over time so that the corresponding hyperparameter is set to zero. Then, we assume the hyperparameters of the source location are same in the three components (east, north and depth). Hence, there are two hyperparameters to be estimated - one for the source location and one for the volume change. We run a two-dimensional grid search and obtain the likelihood contour (Figure 2.3a). The preferred model is given at around (1e1.5,1e0.5). Another take away from the likelihood contour is that a moving source may be not 19 necessary because the likelihood stays almost the same even with a very small location hyperparameter (e.g. 0.1). Figure 2.3. Likelihood (-2L) contours as function of the hyperparameters for volume change and source location (a) when using GPS data only and (b) when using InSAR data only. Star markers represent the preferred model by maximum likelihood estimation. Figure 2.4 shows the estimates of the best fitting model. We can see only small variations less than 200 meters in all the components of the source location. This is consistent with the maximum likelihood results that a relative stable source is indicated. On average, the source is estimated to be in the center of the caldera at ~2.6 km below the sea level. It is notable that the axis ratio is estimated to be ~0.5, which suggest an oblate spheroid source even taking into account the uncertainty. This geometry is quite different from a Mogi source and represents a substantial change in comparison to the pre-eruptive deformation pattern. The estimates of volume change match the temporal features in the time series at OKCE and OKNC. Similar to the pre-eruptive 20 model, the GPS data indicates that the post-eruptive model is spatially stable and the uplifts are mainly due to the volume increase. 50 45 East (m) 0 40 -50 35 Volume Change (106 m 3 ) -100 30 2009 2010 2011 2012 2013 2014 2015 2016 25 20 100 North (m) 50 15 0 10 -50 -100 5 2009 2010 2011 2012 2013 2014 2015 2016 0 2500 2009 2010 2011 2012 2013 2014 2015 2016 Depth (m) 2550 Axis Ratio 2600 0.52 0.51 2650 0.5 2700 0.49 0.48 2009 2010 2011 2012 2013 2014 2015 2016 2009 2010 2011 2012 2013 2014 2015 2016 Time (yr) Time (yr) Figure 2.4. Time-dependent source model based on the GPS data only. All parameters are allowed to vary with time except for the ratio. The solid lines are the estimated values and the error bars show 2σ uncertainties. The reference point of the Cartesian coordinate is the geometric center of the caldera. The volume change is assumed to be zero in the beginning. 2.5.2. InSAR only We invert only a Mogi source for the InSAR data. In other words, we fix the axis ratio to be 1. We first tried to model a simplified spheroid source as we did in the GPS-only inversion. The experiments showed that different initializations return similar estimates and overall goodness of fit and the estimated axis ratio is always close to the initial value. This indicates that the InSAR data have poor constraints on the axis ratio because of lack of near field observations (i.e. observations within the caldera region right above the volcanic sources). 21 The same method as for the GPS data is used to determine the optimal hyperparameters. According to the likelihood contour (Figure 2.3b), the preferred model is given at around (1e2.2,1e0.5). Unlike with the GPS data, the InSAR data suggest a moving source because we do obtain a worse model with a very small hyperparameter for source location. Figure 2.5 shows the estimates of the best fitting model. The model shows that the source moved to the south and to the east by about 1.5 km. This explains why the maximum likelihood solution favors a moving source. Similar horizontal motions were found for this dataset by (Qu et al., 2015) using a series of least squares inversions. In addition, the model shows a slow shallowing of the source. 45 200 40 East (m) -100 35 Volume Change (106 m 3 ) -400 -700 30 2009 2010 2011 2012 2013 2014 2015 2016 25 20 1400 North (m) 15 900 10 400 5 -100 2009 2010 2011 2012 2013 2014 2015 2016 0 2500 2009 2010 2011 2012 2013 2014 2015 2016 Depth (m) 2900 Axis Ratio 1 3300 0.8 0.6 3700 0.4 2009 2010 2011 2012 2013 2014 2015 2016 2009 2010 2011 2012 2013 2014 2015 2016 Time (yr) Time (yr) Figure 2.5. Time-dependent source model based on the InSAR data only. The axis ratio is fixed to be 1. The gray lines are the GPS estimates shown in Figure 2.4. The uncertainties of the GPS estimates are small and not shown here. 22 2.5.3. GPS and InSAR It is obvious that the best fit models for the GPS and InSAR data are not reconciled when inverting for a single volcanic source. Furthermore, we conduct some experiments showing that one spheroid source cannot fit the GPS and InSAR at the same time because there are always jumps in the estimates at the time of the InSAR data acquisitions. We can smooth down these jumps by using smaller hyperparameters but the goodness of fit will be dramatically lower. Recall the results of single inversions, the InSAR data suggests a Mogi source at 3~4 km depth. The source at such depth is consistent with the source determined for the pre-eruptive period, and it was commonly explained as a consistently replenished magma reservoir by previous studies (Lu & Dzurisin, 2014), the post-eruptive source model should also consist of this structure. In Figure 2.5, we compare the InSAR model with the GPS model. It is noticeable that the InSAR model shows a deeper source than the GPS model, but with slightly smaller volume changes. As a result, the InSAR model will most likely underpredict the GPS data, especially for OKCE. We therefore speculate that there should be shallower sources to account for additional and more localized deformation. Therefore, we propose a more complicated source geometry which consists of a Mogi source and a sill source. We use a circular sill to represent a thin layer of magma here because this geometry has the fewest unknown parameters, but it could be a more complicated shape. Even though, it is still much more difficult to model two sources with the data because of the tradeoffs between the source parameters. We take several measures to decrease the degrees of freedom in this scenario. First, we estimate the sill location to be constant rather than varying in time, which mitigates the trade- off between the Mogi source location and the sill location so that we can better separate the two sources. We think this is reasonable for two reasons: (1) we expect the shallow sill to be shallower 23 than 2 km. At such a shallow depth, we think the sill should be a fairly stable structure in nature and separated from the Mogi source; (2) the GPS data should contain more information from this shallow sill and only small variations in time are observed in the location estimates from GPS (Figure 2.4). Second, we run a series of inversions with the Mogi depth fixed ranging from 2.5 km to 4.0 km. We find that the sill depth is coupled with the Mogi source depth.The best model with minimum RMS value is obtained when the Mogi depth is fixed at ~3.2 km. In the final inversion, we thus fix the Mogi depth to be 3250 m. In terms of hyperparameters, we use same values for the sill as that for the spheroid source in the GPS-only inversion and same values for the Mogi as that for the Mogi source in the InSAR- only inversion. Note that the location hyperparameter only applies to the horizontal location of the Mogi source since we fix the Mogi depth. To assign the weights for each dataset, we use the HVCE method described in section 3.4. Because the time intervals of the GPS and InSAR measurements are so different (weekly vs yearly), here we rescale the data covariance based on the residuals of the complete filter run in time and apply a single scaling factor to each dataset. The GPS scale factor is 1.12, indicating that the data noise covariance for GPS is closely true. The scale factors for the T93, T116, E115 and E222 covariance are 1.89, 1.45, 1.77 and 2.12 respectively. Figure 2.6 shows the estimated parameters of the sill source and the Mogi source. The estimates of the Mogi source are similar to the InSAR model. The estimated sill is located to be centered at 53.422 N and 168.123 W, southwest of the geometric center of the caldera, and at ~0.9 km below sea level. The GPS residuals indicates an overall adequate fit of the source model (Figure 2.9a). Considering the InSAR residuals on the volcano flanks outside the caldera region, the model also provides a pretty good fit to the InSAR data (Figure 2.9b). We still calculate the InSAR 24 residuals inside the caldera although we do not use these data. We can see that very large residuals stand out inside the caldera. We believe this is not overfitting but actual deformation associated with the new deposits ejected during the 2008 eruption. More discussions about it are in section 2.7.1. Figure 2.6. Time-dependent source model based on the GPS and InSAR data. (a) Estimates of the deep Mogi source. All the parameters are allowed to vary with time except for the depth. (b) Estimates of the shallow sill. All the parameters are estimated as constant except for the volume change. 25 Figure 2.7. Deformation caused by the Mogi source at (a) GPS sites, (c) TS93 and (e) T116. Deformation caused by the sill source at (b) GPS sites, (d) TS93 and (f) T116. Horizontal and vertical GPS displacements are shown by black and red vectors, respectively. InSAR displacements are shown in line-of-sight direction. This model fit both GPS and InSAR data very well and also explains why the GPS and InSAR data estimate two distinct sources: different coverages in space. The deformation caused by the shallow sill is mainly restricted to inside the caldera, while the deep Mogi source produces a long wavelength deformation pattern with comparable deformation inside and outside the caldera (Figure 2.7). The total deformation is the superposition of these two patterns. All the GPS sites are located inside the caldera or close to the caldera rim (Figure 2.1). It makes sense that the GPS-only inversion returns an oblate spheroid source whose location and shape is in between a sill and a Mogi source. The simulated experiments also support that the combination of a sill and a Mogi source will produce similar deformation at where the GPS sits are. In contrast, the majority of the 26 InSAR coherent pixels lie on the volcano flanks and very few pixels are inside the caldera (Figure 2.2). Thus, the InSAR data do not contain much information from the shallow sill source. It is therefore not surprising that a single Mogi source is sufficient enough to fit the InSAR data. 2.6. Volume change from 2008 to 2019 Another goal of this study is to track the volume change history since after the 2008 eruption. Having the source geometry determined, we then include the GPS data through 2019. From 2016 to 2019, the GPS data are from four continuous sites (OKNC, OKCE, OKSO and OKFG). We fix the source location determined in section 2.3 (Figure 2.11) and estimate the total volume change from both sources. Figure 2.8 shows the estimates of the total volume change. The kinematic volume changes reveal five consecutive exponential inflation episodes from 2008 to 2019. Each inflation episode begins with a rapid increase in volume and then slows down. We fit each of inflation episode using an exponential decay, that is ΔV · (1 − exp (−𝑡/𝜏)), where the characteristic time scale τ and characteristic strength ΔV are the two parameters estimated for each episode. From the first to the fifth inflation event, the characteristic time scale decreases from 1.5 to 0.3 year. From 2008-2019, the total volume change is about 0.068 km3, of which 0.056 km3 is from the deep Mogi source and 0.012 km3 from the shallow sill. The volume change ratio between the deep Mogi source and the shallow sill source is about 5:1 through time. 27 70 Total Volume Change (10e6 m ) 3 60 =0.3 5 V =8 5 50 4 =0.4 V =13 40 2 =0.7 3 V =6 =1.3 4 2 V =19 V =24 3 30 1 1 =1.5 20 10 0 2009 2010 2011 2012 2013 2014 2015 2016 2017 2018 2019 Time (yr) Figure 2.8. Total volume change from the sill and the Mogi source. The exponential fit and 95% confidence are denoted by solid and dashed lines, respectively. The time scale τ and characteristic strength ΔV of each inflation event is shown in the figure. 2.7. Discussion 2.7.1. Subsidence within the caldera region As shown in Figure 2.9b, the InSAR residuals within the caldera region are abnormally large. To compare InSAR with GPS, we project three-dimensional GPS measurements at OKNC and OKCE into the line-of-sight (LOS) directions (Figure 2.10).The GPS data show obvious faster inflation than the surrounding InSAR measurements, especially at OKCE. We believe that the OKNC and OKCE measurements accurately record volcanic deformation as they are both installed on exposed solid rock. But most InSAR pixels inside the caldera are located on ground covered by the 2008 tephra deposits. It appears that these pixels measured not only the volcanic deformation 28 but also significant localized deformation related to the 2008 tephra deposits and loose materials beneath that. Figure 2.9. (a) Residuals of the GPS time series. (b) Residuals of selected interferograms in Figure 2. The linear phase ramp of each interferogram has been removed. The estimated initial position and the common mode has been removed for each time series. 29 Figure 2.10. (a) OKCE monument viewed from Cone E looking NE toward the gates of the caldera and Cone D. The 1997 lava flow is visible on the left. (b) OKNC monument viewed looking NNE toward to Cone B which is visible to the right. The monument is installed in rock, which is covered by ash from the 2008 deposit. (c-f) GPS timeseries projected into LOS direction in comparison with the InSAR LOS timeseries. Random walks at the GPS sites and linear phase ramps have been estimated and removed. One of the possible causes is the self-compaction and erosion of new tephra deposits produced by the 2008 eruption. The caldera and eastern flank was blanketed by thick fine-grained tephra deposits. The thickest part is near the new cone exceeding 25-50 m thickness and thinning to centimeters thickness to the east (J. F. Larsen et al., 2015). The unconsolidated tephra deposits are subject to rapid compaction and are vulnerable to rapid remobilization by strong Aleutian rains and runoff from snow melt; significant erosion was observed in the field in the first few years after the eruption. From September 2008 to August 2010, the accumulation of sediments had enlarged 30 the area of both intra-caldera lakes (J. F. Larsen et al., 2015). Another possible cause is continued viscoelastic relaxation and thermo-elastic cooling from older flows that underlie the 2008 deposits. Localized subsidence inside the caldera was also observed between the 1997 and 2008 eruptions by (Lu et al., 2005) in the region where the historical flows were emplaced. They explained the subsidence as the viscoelastic response of pre-1997 flows due to the loading of 1997 lava flow and the thermo-elastic contraction of 1958 lava flow. These two processes can cause subsidence of centimeters per year in the long term, and there is no reason that this deformation would have stopped after the 2008 eruption. We infer that the InSAR residuals also include the viscoelastic relaxation of the pre-2008 lava flows, especially the 1997 lava, due to loading of the newly emplaced 2008 tephra deposits and the thermo-elastic contraction of the 1997 lava flow. In summary, we think that the large residuals within the caldera in InSAR data are likely sustained non-volcanic deformation, such as self-compaction and viscoelastic relaxation. This is also the main reason why we masked out the InSAR pixels inside the caldera when doing the inversions. 2.7.2. Shallow magma sources A long-lived spherical source at 3 km depth has been widely identified beneath Okmok volcano. It is likely that the Mogi source obtained in this study represents the same source that was found before the 2008 eruption, because of their close depths (Figure 2.9). The shallow sill obtained here is a new structure which was not seen in any pre-eruptive source models. To further prove that this is a substantial difference, we apply our method to the 1997-2008 data set and estimate a spheroidal source model instead of the Mogi source which is usually assumed in previous studies. The best fit axis ratio is 1.1 which is very close to 1. This confirms the findings of the earlier studies that a single Mogi source explained the pre-eruptive deformation. 31 Figure 2.11. Geographic and depth map of the sources. The red star represents the Mogi source and the red dashed circle (line) represents the sill source obtained by combining the GPS and InSAR data. The blue circle is the spheroid source given by the GPS data. For comparison, the black circle denotes the pre-eruptive source. The Mogi source determined by the InSAR data only (section 2.5.2) is close to the Mogi source in the joint inversion (red star). To test whether the differences of GPS dataset between the inter-eruptive and post-eruptive period have an influence on the inversion results, we use a subset of the pre-eruptive GPS data used in (Fournier et al., 2009). They are site OK13, OK31, OKSO, OKCE, OKCD and OKFG. The network made by these six sites has very similar distribution as the network for the post- eruptive period. We see the estimates using the subset of the data is almost identical to using the full dataset but have little larger uncertainties. This proves that even though the sparsity of GPS data for the post-eruptive period, the estimates should not be biased or incorrect. We use synthetic InSAR data to test whether the absence of InSAR observations inside the caldera for the post- eruptive period have an influence on the inversion results. We presume a Mogi source at 3 km 32 depth, which is the estimated depth of the deep Mogi source is in the manuscript, centered in the origin of the x-y coordinate with its volume change increasing exponentially over time. The synthetic InSAR data are generated using the forward model of a Mogi source, a linear phase ramp and random spatial correlated noise. The incidence and azimuth are set to 10o and 340o. We apply our method to the synthetic data then to the data with the observations inside the caldera masked out. We manage to recover the time-dependent source model in both cases but the uncertainties are larger when using the masked data. This proves that we can still obtain the correct estimates even without observations inside the caldera, and confirms that the change in the shallow system found in our study is not because of differences in the spatial sampling between the pre-eruptive and post-eruptive datasets. Note that deflations of sill structures were inferred during the 2008 eruption. Using InSAR measurements, (Lu & Dzurisin, 2010) found closing of a rectangular sill at ~0.5 km depth during the 2008 eruption. They explained this structure as the extraction of groundwater related to the phreatic nature of the 2008 eruption. The shallow sill we obtained might be the same structure given the similarity in depths, but our results indicate that it is more likely a magma body or crystal mush that was drained in 2008 and has been refilling after the 2008 eruption. Freymueller and Kaufman (2010) proposed that some of the 2008 magma had been resident at a shallower depth since before 1997, likely from an older pre-existing magma body that was more chemically evolved (Larsen et al., 2013). The shallow sill source could represent the existing magma body these authors mentioned, which was drained during the 2008 eruption. During the 1997 eruption, Mann et al. (2002)also determined a collapsing sill at ~1.8 km depth along with a Mogi source for the co-eruptive InSAR measurements. 33 We infer that Okmok’s shallow plumbing system is composed of a stable magma chamber at about 3 km depth and some smaller magma storages at shallower depths. This is more in agreement with the petrological constraints (Larsen et al., 2013). Additionally, we infer that some of the shallow bodies are not fed during the inter-eruptive period, so are invisible from geodetic measurements, until they are somehow connected to other sources (e.g. due to eruptions). Our model supports that a shallow thin structure has become a part of the plumbing system and fed since the 2008 eruption. Another scenario is that the shallow magma bodies may not be sustained by the low pressure in the magma chamber after an eruption. For example, the sill found and determined by Mann et al. (2002) during the 1997 eruption might fully collapse as a single Mogi source was still able to well fit the data during the 1997-2008 period no matter what data type and spatial coverage. 2.7.3. Magma supply dynamics An exponential inflation can be explained by a simple hydraulic model made up by a single shallow magmatic source fed by an infinite deep magma reservoir. If the deep magma reservoir is at a constant pressure, the shallow source will inflate exponentially until the pressure difference between the shallow source and deep reservoir vanishes (Pinel et al., 2010). Based on this hydraulic model, our kinematic volume change indicates that 1) the pressure of the deep source may increase episodically and 2) the influx from deep to shallow sources can increase with time at Okmok volcano. First, the pressure of deep source must not be constant in time, which breaks the pressure balance between the shallow and deep sources. More specifically, we interpret that the pressure of the deep source increases in episodic pulses so that each new inflation episode occurs as a result of episodic magma supply into this deep source. Magma supply has been shown to be episodic not 34 only at Okmok but also at other volcanoes. For example, before its 2010 eruption, repeated episodes of magma intrusion and accumulation were found beneath Eyjafjallajökull, Iceland (Sigmundsson et al., 2010). In the Aleutian arc, multiple episodic inflation events have been observed at Akutan volcano from 2005 to 2017 (DeGrandpre et al., 2017; Ji et al., 2017). Episodic inflation events are common at Kilauea volcano in Hawaii as well (Cervelli & Miklius, 2003); in that case, episodic deflation-inflation events were interpreted to result from short-lived blockages in the magma supply. Second, the influx from deep to shallow sources must not be constant in time. Otherwise, all the exponential inflating episodes should have the same time scale. As the exponential relaxation time scale (τ) is decreasing from earlier to later episodes, the influx should somehow increase in time. According to (Lengliné et al., 2008), the influx or time scale is controlled by the radius of the shallow source (𝑎𝑟 ) and width and depth of the conduit connected to the deep reservoir (𝑎𝑐 and 𝐻𝑐 ). Particularly, τ is proportional to 𝑎𝑟3 and inversely proportional to 𝑎𝑐4 . This means that a small change of these two parameters may change the time scale in a significant way. In practice, we do not think 𝑎𝑟 tends to be smaller in time. It is more likely that 𝑎𝑐 increase in time, which would lead to more magma being able to flow through the conduit per unit time and thus a shorter time scale, but other parameter changes are possible. Alternative scenarios include the development and clearing of blockages in the magma transport system (Cervelli & Miklius, 2003), or the formation of additional conduits or dikes for each new pulse of magma from depth. Our present geodetic data are not sensitive to these parameters which limits our ability to test these alternate hypotheses. 35 Figure 2.12. Conceptual model of the magma plumbing system at Okmok and dynamics of deep magma sources. The shallow sill source is inferred to have provided some magma (likely the more evolved component) for the 2008 eruption (Freymueller and Kaufman, 2010; Larsen et al., 2013) but had not been connected to the main shallow magma body prior to the 2008 eruption. Based on a hydraulic model (Lengliné et al., 2008), our results indicate that the pressure of the deep source may be episodic and the magma influx from the deep to shallow source may increase over time. 2.7.4. How close is Okmok to the next eruption? In the current cycle, the volume withdrawal during the 2008 eruption is estimated to be 0.14 km3 based on the InSAR measurements (Lu & Dzurisin, 2010) and 0.11 km3 based on the GPS measurements (Freymueller & Kaufman, 2010) (the value was reported as 0.35 km3 in the original paper but the scaling between the Mogi source strength and volume change was accidentally applied twice; we reran the inversion using the same data and correct the value here). By the beginning of 2019, we estimate that 0.068 km3 of magma has been added since the 2008 36 eruption, which is 50-60% of the volume withdrawn during the eruption. Considering the previous eruptive cycle, the inter-eruptive volume increase from 1997 to 2008 was 85–100% of the 1997 eruption volume (Lu et al., 2010). This may imply the capacity of Okmok magma storage system has a specific range, and an eruption will likely occur when the capacity is reached. Additionally, the average inflation rate declined in time before the 1997 eruption and some oscillations (including episodes of deflation) were observed during 2004 to 2008 (Fournier et al., 2009; Lu et al., 2010). We have not yet seen deflationary features in the current cycle. It is unclear whether these features can be regarded as potential precursors. Future data should help uncover this question. However, kinematic models generally have a limited ability to forecast stability or instability of the magma system. More insights may be obtained by the use of physics-based models (Anderson & Segall, 2011), although those models have more degrees of freedom because of more unknown parameters and requires caution in the basic understanding of the source geometry before using them. 2.8. Conclusion In this paper, we introduce a time-dependent inversion filter combining InSAR with GPS for volcano deformation modeling. We model the post-eruptive deformation at Okmok volcano, Alaska, based on the GPS and InSAR time series observations. We find that a Mogi source together with a shallow sill are required to explain both the GPS and InSAR data at the same time. The Mogi source determined here is close to the general pre-eruptive Mogi source which should represent the long-lived magma body that already existed since before the 1997 eruption. The shallow sill is likely an intrusion into a shallow space that was opened when an older, more evolved magma body was partially drained by the 2008 eruption. We infer that the shallow plumbing system of Okmok is more likely composed of multiple source at different depths. From the end of 37 the 2008 eruption to 2019, the total volume increase is at least 50% of the 2008 volume decrease. We have identified five inflation episodes, each of which has the form of an exponential decay. The characteristic time scale of these inflation events is decreasing with time. We have not yet observed any oscillatory features such as those observed before the 2008 eruption. These features of the volume change curve may help us to understand more about eruptive cycle and unrest at Okmok volcano. We demonstrate that combining GPS and InSAR is crucial to understand the shallow magma plumbing system of volcanoes. In the case of Okmok, the GPS data and the InSAR data individually suggest two distinct sources. We would not be able to figure out the “Mogi+sill” model without combining them. Besides, the InSAR time series techniques are getting more and more popular in volcano geodetic studies, and the time interval of InSAR acquisitions is becoming shorter and shorter with more InSAR sensors collecting data. There are high demands of an efficient workflow to take full advantage of time series observations from both GPS and InSAR. This paper presents a method for combing GPS and InSAR measurements to invert for the kinematics of volcanic source models. As the basis of the method is a type of Kalman filter, it also provides a means of a near real-time tool for volcano hazards monitoring. 38 Chapter 3: 25-year history of volcano magma supply in the east-central Aleutian arc, Alaska 3.1. Abstract Using ~25 years of GPS and InSAR observations and the time-dependent inversion filter, we invert for the volume change history of three active volcanoes in the east-central Alaska Aleutian arc, and the secular velocity of the region. The inferred time-dependent volume change shows 1) two exponentially decaying pulses followed by a nearly linear inflation period at Westdahl volcano; 2) multiple episodic inflation events at Akutan volcano; 3) transient volcanic inflation and long-term deflation signals at Makushin volcano. Comparing our regional velocity estimates with recently published block models, we find a small signal that can be explained by slip deficit on the megathrust, either from a shallow (<20 km depth) almost fully-locked zone or a weakly locked zone confined to 25-50 km depth. 3.2. Introduction The east-central Aleutian arc, Alaska, has several large, active volcanoes, including Westdahl, Akutan and Makushin (Figure 3.1). All three volcanoes formed around the early Pleistocene (Miller et al., 1998). Westdahl is a young glacier-capped shield volcano at the southwest end of Unimak Island. Its edifice is 18 km in diameter, and is composed of a thick sequence of pre-glacial basalt lava flows. Makushin volcano is located on Unalaska island west of Unalaska/Dutch Harbor, the largest city in the Aleutian arc. It is a broad stratovolcano about 1800 m high and 16 km in basal diameter. Akutan volcano is located between Westdahl and Makushin. It is a composite stratovolcano with a circular summit area about 2 km across located in the west- central part of Akutan island. The most recent eruptions of all three volcanoes happened in the early 1990s, producing mostly basaltic lava and volcaniclastic deposits (Miller et al., 1998). 39 Figure 3.1. Tectonic setting of Aleutian arc, Alaska (inset) and topographic map of the study region (black rectangle in inset map) which includes Westdahl, Akutan and Makushin volcanoes. The red lines show the rigid blocks for Alaska Peninsula and eastern Aleutians from Elliott and Freymueller (2020). The brown polygon outlines the 1957 earthquake aftershock zone (House et al., 1981). The focal mechanism represents the 2003 M6.6 Unimak Island earthquake, which happened 200 km east of Dutch harbor just outside the study region. The triangles show the locations of the continuous (red) and campaign (blue) GPS sites. The red vectors denote the average GPS velocity at the continuous sites, with the estimated regional velocity removed. 40 Surface deformation accompanies the volcanic activities at all three volcanoes. Past studies have used measurements of surface deformation (mainly by GPS and InSAR) to study the locations, geometry and temporal evolution of their magma storage systems. At Westdahl volcano, Mann & Freymueller (2003) found that a Mogi source at about 7 km depth could account for the inflation from 1998 to 2001 recorded by a campaign GPS dataset. Lu et al. (2003) used InSAR to study the period 1992-2000, finding that the rate of inflation seemed to decay exponentially with time. Similarly, Gong et al. (2015) processed SAR data for 2006-2015 and also found that the time series of volume change showed an exponential decay. For Makushin volcano, Lu et al. (2002) processed several InSAR images from 1992 to 2000 and found pre-eruptive inflation before the 1995 eruption. However, each of these studies focused only on a specific time window, and did not combine GPS and InSAR data from different satellite missions or viewing geometries. Here, we will use most of the available InSAR and GPS data to recover the magma supply history of the three volcanoes over the 20-25 years since the first geodetic measurements were made there. Besides volcanic activity, the GPS positions include the effects of block motion and subduction. The Aleutian islands are moving in a trench parallel direction relative to the North American plate, toward the southwest, on a microplate called the Peninsula block (Li et al., 2016). The slip deficit on the subduction megathrust varies considerably along strike direction, with a wide region of high slip deficit observed along the eastern part of Alaska peninsula (Li & Freymueller, 2018), transitioning to a creep-dominated megathrust in the eastern Aleutians (Cross & Freymueller, 2008). Given the long time span of the GPS data collection in the region, we can make a more accurate estimation of the secular velocity for the eastern Aleutians than before, if we model it together with the GPS data from the volcanoes. Furthermore, we can evaluate the block models and the subduction coupling models. 41 3.3. Data We use GPS data from 78 sites listed in Table 3.2. Among them, 36 are continuous sites of which 8 were established by the Alaska Volcano Observatory (AVO) starting as early as 2002, and the majority were established as part of the EarthScope Plate Boundary Observatory (PBO). The other 42 sites are campaign sites that were surveyed episodically by the University of Alaska Fairbanks/AVO, starting in 1997. The first campaign measurements in the Dutch Harbor area were made in 1996 by another group (Avé Lallemant & Oldow, 2000). Daily position solutions for the GPS data were estimated in the ITRF2008 frame following the processing procedures described by Fu & Freymueller (2012). We use the time series of weekly averaged position measurements (by GPS week) for the purpose of reducing the computational time. The Westdahl data in the winter were discarded because of large position excursions due to snow and ice accumulation on the GPS antennas at several Westdahl sites. Winter data from the other volcanoes were not affected by this problem. On February 18, 2003, a magnitude Mw6.6 earthquake occurred about 200 km ESE of Akutan island. Based on the USGS focal mechanism (https://earthquake.usgs.gov/earthquakes/ eventpage/usp000br01/moment-tensor), the earthquake ruptured a shallow part of the Alaska subduction thrust. It should produce ENE-directed displacements across Akutan island as a thrust earthquake. There are obvious offsets in the timeseries of sites AKLV and AKMO when the earthquake happened. No other continuous sites had measurements right before and after the earthquake. We estimated the coseismic displacements at AKLV and AKMO by detrending the timeseries and compared the positions before and after the earthquake. Because the whole of Akutan island is far to the epicenter, the displacements at AKLV and AKMO should be almost same so we averaged their coseismic offsets as the final estimate, that is about 10 mm in the east 42 component and -4 mm in the north component. However, these values are about twice as large as the predicted offsets based on the USGS moment tensor solution. After many experiments of adjusting the values of parameters such as epicentral location, dip, strike, etc., we found that the simplest way to match the observations was to increase the seismic moment by a factor of 2. Shifting the event to slightly closer location and using a slightly smaller scaling for the moment would fit the data equally well, and also predicted the same displacement correction at all the sites of interest. A factor of 2 in moment is equivalent to a magnitude difference of ~0.2, so this correction is reasonable considering the uncertainty in seismic magnitude or moment estimates. We then calculated the predictions at all the GPS site locations using the scaled point source solution and corrected the time series observations for all sites, although the corrections are small at Makushin and Westdahl. At Westdahl, the corrections make no detectable difference in our inversion results, while for Makushin the removal of the earthquake displacements reduced the estimated magnitude of the 2001-2004 inflation event at Makushin. Lu et al. (2003) and Gong et al. (2015) processed the InSAR images/timeseries using ERS1/2 and ENVISAT data for Westdahl volcano throughout the period of 1990s and 2000s. For Makushin volcano, Lu et al. (2002) and Lu & Dzurisin (2014) processed the InSAR images using ERS1/2 and ENVISAT data. We directly use these published InSAR data for our study in the form of InSAR time series. Readers are referred to the original papers for details about data processing. We do not use any InSAR data for Akutan volcano because the GPS network at Akutan has consisted of 12 continuous sites distributed well across the whole island since 2005, which means that the GPS data alone are able to provide good constraints on the deformation source parameters, and long-term InSAR coherence was relatively poor. 43 3.4. Method Previous studies have shown that a spherical “Mogi” source is sufficient to model the volcanic deformation at Westdahl and Makushin volcanoes (Gong et al., 2015; Lu et al., 2002, 2003; Mann & Freymueller, 2003). Akutan has a more complex magma plumbing system (DeGrandpre et al., 2017; Wang et al., 2018). DeGrandpre et al. (2017) observed an asymmetry in the displacement pattern at Akutan over certain time intervals, but they modeled velocities relative to a specific PBO site, and implicitly assumed that the regional block velocity was equal to the velocity of that site. That assumption does not match other observations from nearby islands, but some asymmetry cannot be ruled out. They used spheroid/ellipsoid sources and found that the best models, in terms of source geometry, for different time periods were different. Ji et al. (2017) extracted the displacements during four major inflation events from 2005 to 2017 and used a Mogi source to model those displacements. They removed a long-term trend from the data as the first step of their analysis, so they modeled only a portion of the transient events (their de-trending also removed some volcanic signal, as shown by DeGrandpre et al. (2017)). We tried to fit the Akutan data using a spheroidal source and then compared the results with those using a Mogi source. We used the same hyperparameter that controls the temporal smoothness of the volume change for both cases. The overall reduced Chi-squared is 1.97 when using a Mogi source and 1.91 when using a spheroid source, and the F test (Dzurisin et al., 2009) does not support the hypothesis that a spheroid source provides a significantly better fit than a Mogi source. According to the results from past studies, we think that the major inflation events at Akutan volcano can be explained well by the overpressure of a Mogi-like source, while between these major events different source geometries might fit the data better. One interpretation could be that different parts of the magma system are activated during the major events and between them. However, the volume change 44 caused by the activity between the major inflation events are very subtle and smaller in amplitude than the inflation during events. Therefore, our model geometry estimate is dominated by the major inflation events and we do not consider the complexity of the source geometry at Akutan for the purpose of studying the volume change history here. In summary, we found a Mogi source is still sufficient enough for studying the volume change history of Akutan volcano because the largest inflation events were well described by a Mogi source. Therefore, we assume Mogi sources that are not moving over time (i.e. the source location is estimated as a constant in time) to represent the pressurized magma chamber producing the surface deformation for all three volcanoes. Besides the location, the Mogi model features a source strength parameter that combines volume and pressure change (inseparably). Delaney and McTigue (1994) showed that this parameter can be related to the injection volume, which they defined as the change in the volume of the walls of the magma reservoir, and this is the measure of volume that we use for the rest of this paper. In the case of incompressible magma, this would also be the volume (proportional to mass) of magma intruded into the shallow storage system from greater depth, but we have no tight constraints on the magma compressibility at these volcanoes, so we cannot convert injection volume to magma mass change without additional assumptions. To invert the time series, we use a time-dependent non-linear inversion filter based on the unscented Kalman filter (Xue et al., 2020). At each time step of the inversion filter, the source parameters and their covariance (i.e. volcanic source locations and volume changes, and the regional secular velocity relative to North America) are updated based on the new data and a state transition matrix. The source location is estimated as constant in time, while the volume change for each volcano is allowed to vary with time. The reference volume change is based on the earliest measurement (GPS position or SAR image) for each volcano. The regional secular velocity 45 parameter is assumed to apply to all volcanoes, and only the GPS data are sensitive to this parameter. Other additional nuisance parameters are estimated for each data set as described in section 2.4. A backward smoother is used to smooth the estimation after forward filtering. 3.5. Results The estimated source parameters, including location and time-dependent volume change, for all three volcanoes are shown in Figure 3.2. Westdahl volcano shows two distinct exponential inflation episodes followed by linear inflation (Figure 3.2a). The first exponential inflation started after the 1992 eruption and ended by 2001, when the second exponential inflation initiated. Starting from 2007, the volume has increased almost linearly with time. Previous studies analyzed shorter periods of time, and we can compare our results with the volume changes inferred over those periods of time (colored lines drawn in Figure 3.2a). Our estimates provide more temporal detail and continuity, but are overall consistent with the previous studies during each study’s period of data. In our results we can clearly see the sudden change in the volume change rate in 2001 due to the initiation of the second exponential inflation event, which was not clearly identified in any of the previous studies as it started in the gap between the first two InSAR data sets (ERS-1/2 to Envisat). Campaign GPS measurements across this time period help to constrain the onset of this inflation event, but only when combined with the full time series of all data (including 2014-15 surveys of a few campaign sites). Mann & Freymueller (2003) only had a short time span (3 years, 1998-2001) of the campaign GPS data available, and while they correctly found an average inflation rate that was faster than the late 1990s InSAR rate (Figure 3.2a), that data by itself could not constrain the onset of the second exponential inflation event. 46 Figure 3.2. Time-dependent volume change estimates with 2σ uncertainties (error bars) for (a) Westdahl, (b) Akutan and (c) Makushin. The gray curves or lines are fit to the estimates. Other colored lines show the estimates from previous studies. The caption for each subplot also lists the 47 Figure 3.2 (cont’d) estimates of the location (longitude, latitude and depth in km) of the Mogi source for the corresponding volcano. The uncertainties for the source location estimates are ~100 m in horizontal and ~200 m in depth, except the depth for the Makushin source which is fixed to be 7 km. Akutan volcano has shown multiple episodic inflation events that have occurred approximately every 2-3 years (Figure 3.2b). The inflation events in 2008, 2011, 2014 and 2016 already were identified by previous studies (DeGrandpre et al., 2017; Ji et al., 2017; Ji & Herring, 2011), and the estimated volume increase of each event is consistent with the previous results. We also find two earlier events that happened in 2004 and 2006 when only the AVO continuous GPS sites were operational and an inflation event that has been ongoing since 2019. Each of these episodic inflation events can be fit by an exponential decay time history, and their time scales range from 0.3-0.8 yr. Before 2002, only campaign data were available, but one or more inflation events must have occurred between 1997 and 2002, although only the total inflation over that time can be determined. The average rate of inflation over 1997-2002 was similar to that over the later time period. At Makushin volcano, the rapid inflation before the small explosive eruption in 1995 (Lu et al., 2002) and a period of constant deflation from 2004 and 2009 (Lu & Dzurisin, 2014) are both reproduced here (Figure 3.2c). In addition, our estimates indicate a period of deflation from 1996 to 2001 after the eruption. There is another inflation period from 2016-2018, and its volume increase is similar to the 2001-2004 inflation but significantly smaller than the inflation before the 1995 eruption. 48 Figure 3.3. Volume change estimates for Makushin volcano in the bestfit model and the models when the east component of the regional velocity is fixed to the values deviated 1 or 3 times of the estimation uncertainty from the original estimate in the paper. The two black boxes denote the time windows that are affected by the parameter tradeoff between the regional velocity and volume change because of a lack of InSAR or near-field GPS data. Between the two deflating periods, there is period of inflation throughout 2001-2004 although there is a strong tradeoff between the east component of the regional block velocity and the volume change estimate because of the poor distribution of the GPS sites at that time. Before the three continuous sites (MAPS, MREP and MSWB) were installed at Makushin in 2011 and 2012, most of the GPS data were from campaign sites that were all located around Dutch Harbor, well east of the volcano deformation source. One campaign site (MHAA), measured in 2003, 2011 and 2012, was located close to MAPS. Thus, there is a strong tradeoff between the estimates of 49 the volume change and the regional block velocity (mostly in the east component). Specifically, it is hard to tell whether an eastward motion trend in the GPS data from Dutch Harbor is due to a volume increase of the magma storage or just part of the regional velocity. However, small time variations can be seen in the east component time series of the longest-surveyed campaign site in Dutch Harbor (DCH1) and in the AV09 and DUTC continuous sites. We assume that these variations reflect changes in the volcanic source rather than tectonic transients, because a tectonic transient would be expected to affect mainly the north component (the direction of the elastic slip deficit signal, see Figure 3.4) while a volcanic signal would appear mainly in the east component. To examine the impacts of this tradeoff and the uncertainties of the inferred regional velocity, we inverted for the time-dependent volume change for Makushin volcano while fixing the regional velocity to several alternate values corresponding to our best estimate of the regional velocity and the ±3 values. Figure 3.3 shows that the volume change estimates during the period of 2001-2004 and 2011-2012 are obviously affected by the regional velocity assumed. This is expected because those are the two time periods when there was a time gap between separate InSAR stacks, and when GPS data were only available from Dutch Harbor in the far field. Outside of those time periods, the volume change estimates do not change substantially even if the regional block velocity is perturbed because there are enough near-field continuous GPS data or InSAR data to constrain the volcanic source parameters and break this parameter tradeoff. The main question we want to answer here is whether the estimated inflation from 2001 to 2004 is realistic, or could it just be an artifact of the data distribution? According to Figure 3.3, while the amount of the inflation is sensitive to the regional velocity estimate, some volcanic inflation is always present in all models. 50 Our secular velocity estimation indicates that the study region is moving in a primarily trench parallel direction with a rate of ~5 mm/yr (Figure 3.4). Compared with two previously published block models (Elliott & Freymueller, 2020; Li & Freymueller, 2018) (Table 3.1), the east component estimates all agree to within 0.2 mm/yr, but there is about a 2 mm/yr difference in the north component between our estimate and the block model estimates (top left white inset in Figure 3). We will discuss this more in section 3.6.2. Table 3.1. Estimated regional velocities relative to North America in comparison with block model estimates. This study Li et al., (2018) Elliott et al., (2020) East (mm/yr) -4.34  0.35 -4.56  0.35 -4.47  0.63 North (mm/yr) -0.52  0.25 -2.02  1.47 -3.23  2.38 3.6. Discussion 3.6.1. Magma Supply Lu et al. (2007) has shown that the Aleutian volcanoes present diverse patterns of spatial deformation. From our time-dependent volume change estimates, we can see the magma supply dynamics are also different among the three volcanoes studied, despite their neighboring locations (50-100 km). This diversity is also found in some other volcanic arc systems. For instance, Henderson & Pritchard (2013) found significant variation in deformation rates from InSAR observations throughout five volcanoes in the Central Andes Volcanic Zone. Similarly in Kamchatka, different plumbing system geometries were inferred from seismic tomography beneath three volcanoes located close together in an area of approximately 50 × 80 km (Koulakov 51 et al., 2017). These findings reflect that the magmatic systems feeding arc volcanoes are generally diverse in time and space. Dzurisin and Poland (2019) stated that the magma supply rate is likely a primary control on the eruptive behavior and hazards of Kilauea volcano. Variations in magma supply rate in Aleutian volcanoes also have been observed in the past, as well as in this study. Lu et al. (2003) used a simple model for the hydraulic connection between shallow and deep sources in the presence of an initial pressure difference (e.g. Pinel et al., 2010) to explain the exponential decay in magma supply rate at Westdahl volcano in the first several years after since the 1992 eruption. In the context of that model, the inflation cycle is completed when the pressure difference has decayed to zero. This model can explain only a single exponential decaying inflation (EDI) in the volume change. However, both Westdahl and Akutan have experienced several EDIs so far (Figure 3.2). Multiple EDIs have also been observed at Okmok volcano, which is about 50 km southeast to Makushin volcano, during its inter-eruptive periods (Fournier et al., 2009; D. Walwer et al., 2019; Xue et al., 2020). Walwer et al. (2019) derived a model considering the viscosity gradient in the vertical conduit and magma supply, that in certain condition can result in temporal variations in the volume change. Their model requires the formation of a new shallower magma body for each EDI. An alternative explanation is that the pressure of deep magma sources increase in a stair- stepping way in time rather than being constant due to the arrival of discrete pulses of new magma (Xue et al., 2020), and we favor this explanation. Besides the EDIs, it is also notable that Westdahl transitioned to nearly linear inflation around 2007, although that could be an EDI with a very long characteristic timescale. No matter which case, our results indicate changes of the magma supply over time. Further work is needed to investigate whether and how these variations/changes relate to eruptive hazards. 52 Another measure to evaluate the volcano hazards is the total volume change recovered from the previous eruption. Given that the total volume increase at Okmok between the 1997 to 2008 eruption was 85-100% to volume loss in the 1997 eruption, Lu et al. (2010) proposed that the timing of next eruption of Okmok might be when the volume increase compensates for the volume loss in the previous eruption. Among the three volcanoes in this study, Westdahl had an eruption in 1992, and Lu et al. (2003) estimated the volume decrease of that eruption to be ~0.05 km3. However, the total volume change since then (~0.16 km3) is over 3 times larger than the co- eruptive volume decrease. The 1995 eruption of Makushin was very small, and much smaller than its pre-eruptive inflation or the volume changes observed since 1995. It is clear that a simple “volume predictable” eruption threshold model does hold for these volcanoes. We cannot evaluate Akutan at this point because there is no documented estimation on the volume loss of the shallow magma storage for its last eruption. There was a small explosive eruption at Makushin in 1995, and it was followed by a period of deflation. Overall, the volume change history of Makushin reveals several inflation-deflation patterns, in contrast with Akutan and Westdahl that are inflation dominated. We think these long- term deflations in the absence of an eruption or intrusion are most likely a volume/pressure decrease related either to degassing or cooling/crystallization (Lu & Dzurisin, 2014), and they are not unique in the Aleutians. Similar long-term deflation has also been observed at Fisher caldera, located near Westdahl, and was explained to due to magma cooling and crystallization as well (Mann & Freymueller, 2003), although that was at much shallower depth. Deflation of the magma body inferred in this study at 7 km depth appears to have occurred after each intrusion, and the inflation and deflation are of similar size. One possibility is that the inflation/deflation in this case is mainly the result of changing gas pressure, with little new magma added. Another is that the 53 deflation is associated with the movement of magma to a much shallower body (<< 7 km) that is not easy to see with the geodetic data. Changes in very shallow magma storage could be mostly invisible to our geodetic data because of the loss of InSAR coherence over the glacier cap and the distance to the GPS stations (~10 km), which mean that there are no data very close to the summit. 3.6.2. Plate Motion and Slip Deficit on the Megathrust As mentioned above, there is a ~2 mm/yr difference in the regional velocity (Figure 3.4) between our estimate and those from two previous block modeling studies (Elliott & Freymueller, 2020; Li & Freymueller, 2018), and the difference is primarily perpendicular to the trench (our estimate shows northward motion relative to the block models). Note that both block models are mainly constrained by data from the Alaska Peninsula region, east of our study area. Therefore, this difference could be due to elastic deformation from the slip deficit on the locked subduction interface. We first consider past slip deficit models (Cross & Freymueller, 2008), which assumed a deep seismogenic zone extending from 25-50 km depth offshore of Dutch Harbor and creeping zone near the trench, corresponding to the aftershock zone outline of the 1957 Mw8.6 Andreanof Island earthquake at the eastern end of the rupture (see Figure 3.1, inset). However, the true extent of the rupture zone for the 1957 earthquake has long been uncertain. Nicolsky et al. (2016) showed that such a deep rupture could not explain the observed high tsunami runup on the Pacific coast of Unalaska island, nor the tsunami arrivals at Dutch Harbor. Therefore, we also consider a model in which there is only a shallow locked zone, starting at the trench and extending from 0-15 km depth, which corresponds to the best-fitting model of Nicolsky et al. (2016). We use a least squares method to estimate the slip deficit rate. We assume that the block model estimates, constrained by nearby data but not including the sites affected by volcanic 54 deformation, contain the regional block velocity only while our data also include elastic strain from the slip deficit. The signal from the slip deficit is small at sites located along the volcanic arc in any models and it is fairly uniform across the study region, so it is absorbed into our block velocity estimate rather than biasing the volcanic source parameters. In this analysis, the observations are velocity estimates from our study, and the block model based estimates of Li & Freymueller (2018) and Elliott & Freymueller (2020), which are shown in Table 3.1. The two block model estimates were based on data from nearby parts of the Peninsula block, so they are independent from our estimated velocity. The key assumption in this inversion is that our estimate contains both the block velocity itself and the effect from the slip deficit, while the two block model estimates do not contain the slip deficit. The unknown parameters to be estimated are the slip deficit rate on the assumed locked zone, and the velocity (east and north components) of the rigid block. The relationship between the unknown parameters and the observations are given in below, 𝑑𝑒1 𝑔𝑒 1 0 𝑑𝑛1 𝑔𝑛 0 1 𝑣𝑠 𝑑𝑒2 0 1 0 2 = 0 𝑣 [ 𝑒] 𝑑𝑛 0 1 𝑣𝑛 (4) 𝑑𝑒3 0 1 0 [ 0 0 1] [𝑑𝑛3 ] where 𝑑𝑒1 and 𝑑𝑛1 are the regional block velocity estimates from this study; 𝑑𝑒2 and 𝑑𝑛2 are the estimates of (Li & Freymueller, 2018); 𝑑𝑒3 and 𝑑𝑛3 are the estimates of (Elliott & Freymueller, 2020); 𝑔𝑒 and 𝑔𝑛 is the Greens function for the velocity due to the slip deficit; 𝑣𝑒 and 𝑣𝑛 are the east and north components of the regional velocity based on the estimates from the three studies; 𝑣𝑠 is the slip deficit rate on the subduction. The coupling coefficient can be then calculated by dividing 𝑣𝑠 by the relative rate of the plate motion (66 mm/yr). 55 In the case of the deeper locked zone, we estimate the slip deficit rate to be 8.1 mm/yr, which represents ~12 % of the convergence rate between the two plates. This value matches the estimate given by Cross & Freymueller (2008). In the case of a shallow locked zone with the deeper part of the interface creeping, the slip deficit rate is estimated to be 61.3 mm/yr, indicating a highly locked zone. We favor the model with only a shallow locked zone, based on the recent slip constraints for the 1957 earthquake determined by Nicolsky et al. (2016). However, because our data are restricted to the volcanic arc far from the trench, and with a narrow trench-normal aperture for the network, the two slip deficit models give equally good fits to the data. New observations from sites closer to the trench, on the Pacific coast of the island or from seafloor geodesy (Bürgmann & Chadwell, 2014), will be required to distinguish between these two models. Figure 3.4. Two alternate plate coupling models: (a) slip deficit only at shallow depth (<15 km) or (b) slip deficit only at 35-50 km. The dashed line outlines the Alaska trench. Slab contours are shown in black lines, and the outline of the locked zone as a red rectangle with the slip deficit rate given as a percentage of plate motion rate. The top left inset denotes the secular velocity from our estimates (red), previous block models (green) and their difference (blue). The percentage number 56 is the slip deficit rate relative to the convergent rate between the Pacific plate and Peninsula block (6.6 cm/yr). The predicted deformation caused by the slip deficit is denoted as black vectors. 3.7. Conclusion Using long-term GPS and InSAR measurements, we estimate the location and the time- dependent volume change during the last ~25 years for the shallow magma storage beneath Makushin, Akutan and Westdahl volcano in east-central Alaska Aleutian arc. The inferred volume change history indicates that the magma supply dynamics of the three volcanoes are diverse: 1) two exponentially decaying pulses followed by a nearly linear inflation period at Westdahl volcano; 2) multiple episodic inflation events at Akutan volcano; 3) transient volcanic inflation and long- term deflation signals at Makushin volcano. We also invert for the secular velocity for the study region and find a small (~2 mm/yr) northward component relative to the velocity estimated from block models. We think this small signal is likely due to elastic strain from slip deficit on the megathrust. Our preferred coupling model is a shallow (0-15 km) nearly fully locked zone, although an alternate model with a low slip deficit rate (~12% of the plate motion) on only a deeper locked zone at 35-50 km depth would also fit the data. 57 Table 3.2. GPS site locations. Name Longitude Latitude Type AC10 -164.887 54.523 continuous AKGG -165.994 54.198 continuous AKHB -165.734 54.104 campaign AKLV -165.956 54.163 continuous AKMO -166.012 54.09 continuous AKPS -165.85 54.145 campaign AKRB -166.071 54.129 continuous AKSO -165.944 54.138 campaign AKU -165.763 54.143 campaign AKUA -165.606 54.144 campaign AKUN -165.593 54.211 campaign AMRG -166.043 54.067 campaign AV06 -165.766 54.147 continuous AV07 -166.039 54.163 continuous AV08 -166.028 54.136 continuous AV10 -165.934 54.098 continuous AV12 -165.898 54.211 continuous AV13 -165.898 54.153 continuous AV14 -165.842 54.119 continuous AV15 -165.71 54.1 continuous AV25 -164.779 54.53 continuous AV26 -164.58 54.572 continuous AV27 -164.723 54.492 continuous AV29 -164.586 54.472 continuous BROD -165.869 54.111 campaign DCH1 -166.528 53.889 continuous DT16 -166.531 53.874 campaign DT19 -166.533 53.874 campaign 58 Table 3.2 (cont’d) DTOC -166.536 53.876 campaign DUTC -166.548 53.905 continuous FLOW -166.078 54.169 campaign FTOP -165.99 54.116 campaign GUNN -166.516 53.924 campaign HSB -165.914 54.186 campaign ILIU -166.485 53.852 campaign LVA2 -166.016 54.164 campaign MAPS -166.94 53.808 continuous MHAA -166.939 53.808 campaign MREP -166.748 53.81 continuous MSWB -166.788 53.915 continuous NFLW -166.035 54.184 campaign OPEN -165.963 54.193 campaign REF2 -166.097 54.114 campaign ROTK -165.521 54.064 campaign SBS2 -166.515 53.897 campaign SCAP -164.745 54.398 campaign UPMO -166.012 54.09 campaign WEFC -164.545 54.613 campaign WESE -164.542 54.538 campaign WESN -164.581 54.571 campaign WESP -164.683 54.558 campaign WESS -164.724 54.479 campaign WFAR -164.779 54.533 campaign WHL -165.826 54.119 campaign WPOG -164.746 54.596 campaign WWLL -166.081 54.15 campaign 59 Chapter 4: Single-station detection of transient deformation in GPS time series based on machine learning methods 4.1. Abstract The development of GPS networks within the past two decades has provided a large volume of dense and continuous time series data, which brings up an opportunity but also a challenge for utilizing the large volume of data. Inspired by recent studies on detection of transient signals in geodetic time series data (e.g. seismic, sea floor pressure) using machine learning methods, we present a simple model for transient detection based on a recurrent neural network structure. Unlike most previous studies that are based on fixed time windows or metrics computed from the time series, our model uses each single data point of the entire time series as input sequentially and directly outputs the transient probability for all points in the time series. The model is trained using synthetic data time series that are relatively short in length, but can be applied to geodetic time series of arbitrary length, and could even be used in a real-time sequential mode. We test and validate the model performance by applying it to a synthetic testing dataset and a real-world dataset for the Cascadia region. For the latter dataset, our detection results are consistent with previous studies on the slow slip events in the Cascadia region. We also test the model on time series containing coseismic offsets and longer-term transients to test how its performance degrades when it encounters data with different characteristics from the training data sets. 4.2. Introduction There has been a tremendous development in automated processing and open data from the Global Positioning System (GPS). For example, the Nevada Geodetic Laboratory currently obtains geodetic GPS data from more than 17,000 stations and updates the daily position coordinates of 60 more than 10,000 stations every week (Blewitt et al., 2018). Such big datasets provide new research opportunities for geophysical problems but also a challenge in terms of how to efficiently extract signals of interested. In this case, we explore the detection of transient deformation due to various geophysical processes. During the past two decades, many methods have been proposed for detecting transient deformation in GPS timeseries. Examples include but are not limited to methods built upon some kind of Kalman filter (Bekaert et al., 2016; Ji et al., 2017), decomposition analysis (Dong et al., 2006; Damian Walwer et al., 2016) or template matching (Riel et al., 2014). These methods usually require a network of GPS stations for the spatial information to get better and more reliable results. To obtain more flexibility, Crowell et al. (2016) proposed a single-station detection method based on the Relative Strength Index (RSI). They used the RSI, an indicator of the divergence from and beyond the norm of the time series, to calculate the transient probability of a data point using a probabilistic approach. This single feature, the RSI value, was demonstrated to be simple but efficient for detecting transient signals in GPS time series. One question would be whether or there are more such metrics like the RSI value that could be used for the task, or if a machine learning (ML) approach can be based directly on the time series values themselves. More recently, machine learning methods have been used as a feature extraction tool for data mining in many geophysical fields, such as phase picking in seismic waves (Ross et al., 2018; Zhu & Beroza, 2018), volcanic deformation from satellite images (Anantrasirichai et al., 2018), and transient detection in sea floor pressure data (He et al., 2020). In the above methods for time series data, the common solution is to scan the time series through a time window and then do the classification just using the data within that time window. This solution does not actually take the 61 advantage of global information in the entire time, series which could be exceptionally useful for classification. In this manuscript, we present a ML model based on the recurrent neural network (RNN) for transient detection in geodetic timeseries in a sequential manner. First, we introduce the structure of our ML model as well as the training process. We then conduct some synthetic tests to evaluate the model performance. Last, we validate our model by applying it to a real-world dataset for the Cascadia region, and then test it against additional data with somewhat different kinds of transient signals that are not appearing during the training process. 4.3. Methods 4.3.1. Synthetic data generation and labeling We use synthetic data to train the model for two reasons. First, a large well-defined labeled database for transient signals in the GPS time series is not accessible at this point. Besides, the transient signals are rare in nature, which limits the total number of the positive samples we could use for training. Second, the noise characteristics and the kinematics of most transient signals are well known, so we can easily generate synthetic data that resembles the real-world data very closely. Therefore, we can obtain as many close to realistic samples as we need for model training using synthetic data. The synthetic daily GPS time series consists of secular, seasonal, transient deformation, and colored noise (Figure 4.1a). The secular signal is a linear trend over time. The seasonal signals are a combination of an annual and a semi-annual sinusoid. For seasonal terms, we use the interquartile ranges of seasonal model terms from the Scripps Orbit and Permanent Array Center (SOPAC) archive for all western North America stations, 1630 stations in total. We consider three types of noise that are commonly found in geodetic timeseries observations - white, flicker, and 62 random walk noise. The amplitudes of each noise component are randomly sampled from the interquartile ranges for Southern California stations (Langbein, 2008). The specific sample range of noise and seasonal terms is the same as listed in the Table 2 of Crowell et al. (2016). Last, we simulate the transient signal using an arctangent function described by three parameters - size, timescale, and center. Figure 4.1. An example of (a) raw synthetic timeseries, (b) detrended synthetic timeseries and (c) the labels for the timeseries. The next step is to label the synthetic data. In general, the label of all data points within the duration of a transient is set to 1, and otherwise 0 (Figure 4.1c). By doing so, the timing of a transient is also determined. Additionally, we set a threshold for the minimal transient size only above which we set the label to be 1. The reason for this is to build a clear decision boundary for 63 the model to learn. Considering that a very small transient is just very close to a pure noisy data in nature, it would be harmful for learning if they are labeled differently. In this study, we label the data point to be 1 only when the signal-to-noise ratio is greater than 1 for the transient signal. In other words, we are targeting to detect transients that have at least the same order of magnitude as the noise level, equaling ~3 mm for the vertical and ~1 mm for the horizontal. 4.3.2. Model structure, training and prediction Figure 4.2 shows the model structure. It first consists of a bi-directional long short-term memory (LSTM) layer with 32 units and a dense layer with 16 units taking each data point in a time series as input chronologically. We use a hyperbolic tangent activation function for the LSTM layer and an ReLU (Rectified Linear Unit) activation function for the dense layer. This results in 9971 parameters in total. LSTM is a special type of RNN that solves the vanishing gradient problem (Hochreiter & Schmidhuber, 1997), thus is the most suitable structure for long time series analysis. The final output layer is a 1-unit dense layer with a sigmoidal activation function. The model can be built more complicated, such as by increasing the units of LSTM layer, and still trainable using more data, but we find that this simple structure already performs well. The total binary cross entropy is used as the loss function, − ∑𝑛𝑖=0[𝑦𝑖 𝑙𝑜𝑔(𝑝𝑖 ) + (1 − 𝑦𝑖 )log (1 − 𝑝𝑖 )], where 𝑦𝑖 and 𝑝𝑖 represent the true label and the predicted probability of each data point from all training time series samples. We adopt the ADAM algorithm (Kingma & Ba, 2015) to minimize the loss. We generated 10,000 timeseries, each with 500 epochs. Each synthetic timeseries has a random number of transients from 0 to 8 (0 means the timeseries just includes noise). The transient timescale ranges from 5 to 30 days, and is sampled uniformly during the data generation. All the transients, if more than one, are randomly spread out without overlapping in terms of their centers. The transient size is uniformly sampled up to 20 mm positive or negative. 64 These parameters can be chosen differently for transient signals with different properties in other cases, but were chosen here to be representative of short-term slow slip events. Figure 4.2. A diagram showing the structure and workflow of the ML model. Each input dn is a data point of a timeseries. The bidirectional LSTM layer and the dense layer are used to extract features. The final output gives the transient probabilities. All the timeseries are then detrended and standard scaled between 0 and 1 (Figure 4.1b). The scaling is helpful to speed up the training process and makes the model independent of the noise level so that it can be applied to all three components of the GPS observations in the same way. We use the same length of synthetic data to simplify the training process, because only cases with the same length of input data can be trained in a single batch. However, the model does not specify the length of the time series so it can be applied to timeseries of any lengths for prediction. By storing the hidden states of the LSTM cells, the model cuold be used for unidirectional prediction for a real-time data stream. 65 The output of the model is the transient probability of the current data point, which is a continuous number from 0 to 1. The higher the probability is, the more likely that the data point is during a transient. A decision threshold is needed to convert continuous probability to discrete labels (either 0 or 1). If the probability is greater than 0.5, the predicted label will be 1, and vice versa. In the synthetic tests, we simply use the threshold of 0.5 to obtain the predicted label from the predicted probability for computing the model performance metrics, although there are more thorough ways to determine the threshold such as by using the ROC curve (Bradley, 1997). For figures that involve model prediction results, we show the direct output from the model, that is the probability. 4.4. Results 4.4.1. Testing on synthetic data We consider two metrics to evaluate the model performance: precision and recall. Recall represents the percentage of all detections that actually are transients, and precision represents the percentage of all transients that are successfully detected. For both metrics, the higher the better. They are defined as follow: TP TP Precision= , Recall= TP+FN TP+FP (5) where TP, FP and FN represent the number of true positive, false positive, and false negative. To illustrate how the model performance changes for detecting transients of different sizes, we conduct a series of tests to obtain the precision and recall score as a function of offset size. We first generate pure noise timeseries as described in section 4.3.1. Next, we add a transient term function of size 1 mm at the 1-year mark in the time series. We then repeat the same procedure but increasing the transient size up to 50 mm in 1 mm increments. A total of 10000 samples are generated to test each size level. Note that in these tests, it will be a TP if there are any detections 66 within the duration of a transient because in reality one detection is enough to correctly remodel the timeseries. It will be a FN if there is no detection at all within the transient. Outside of the transients, all data points are treated independently when counting the number of FP. Figure 4.3. (top) Precision and (bottom) recall curve as a function of transient size. Figure 4.3 shows the score of precision and recall as a function of transient size. The noise characteristics of the north and east components are similar, so we group them together as the horizontal component. From the precision curve, we can see that the horizontal component starts at a precision score of 0.1 when the offset size is 1 mm and rises to above 0.95 by about 3 mm. The precision score of the vertical component takes longer to goes up above 0.95 when the 67 transient size is about 8 mm, which is expected because of the larger noise level. After that, both components have a precision score of above 0.95 with larger offsets. The recall curve looks similar but slightly different. Starting at 1 mm, nearly zero transients are correctly being detected in both components. In the horizontal, the score reaches to 1.0 by about 5 mm and in the vertical by around 15 mm. That is, above those sizes, a detection is guaranteed to occur. Combining with the results from the precision scores, we can conclude that, in the horizontal (vertical), a transient larger than 3 (8) mm of size will be most likely being detected although some false detections will persist, and we are more confident that the detection is correct when the transient size rises up to 5 (15) mm. These testing results show that our ML model could outperform the RSI method which was tested under similar conditions (Crowell et al., 2016) for detecting transients of small sizes. It indicates that more features are extracted by the ML model for make the decision. 4.4.2. Application to a real dataset To validate our ML model in a practical way, we apply it to detecting transient events for the Cascadia region. We use the same dataset that was used in (Crowell et al., 2016) because it has been carefully preprocessed and is easy to obtain from the SOPAC archive (ftp://garner.ucsd.edu/pub/timeseries/measures/ats/WesternNorthAmerica/). The displacement time series was detrended with all model terms including seasonal terms, coseismic offsets, postseismic parameters, and nontectonic offsets primarily due to instrument or antenna changes. Common mode errors were estimated using principal component analysis and removed. Therefore, the data are assumed to contain only transient signals and noise. Readers should refer to (Crowell et al., 2016) for a detailed description of the dataset and their processing workflow. Our model is 68 not trained for the presence of significant volcanic deformation, so we exclude all stations nearby Mt. Helens, leaving a total of 196 stations. Figure 4.4. (top) The east component of displacement and predicted transient probabilities for station ALBH on southern Vancouver Island. (bottom) Transient detections for all stations. Shown are the transient probability greater than 0.5 in any of the three components (east, north, up). 69 We filled data gaps with zeros to ensure a regular time interval but masked out those epochs when predicting. As an example, we present the east component of displacement and the transient probabilities derived by our ML model for station ALBH (Figure 4.4). It is apparent that the spikes in the transient probabilities coincide with the offsets that can been directly seen in the timeseries. We want to point out that there are totally 4015 epochs for station ALBH and our ML model is trained only using synthetic samples contained 500 epochs, which proves that a universal pattern for transient detection has been learned by our model. We followed the same process for all the stations. The transient probabilities that are greater than 0.5 in any of the three components are shown in Figure 4 as a function of latitude and time. It becomes apparent that some detections are coherent in space and time, representing the well known slow slip event in the Cascadia region. The pattern and migration of these events are consistent with those reported by (Michel et al., 2019; Schmidt & Gao, 2010), especially the relatively larger ones in the northern part, such as the two events in 2011. Last, we remodel each input time series using the ML detection results to separate the average rate between slow slip events from the displacements in the slip events. Each time series is remodeled as a linear combination of a constant at time zero, a velocity, and a set of Heaviside step functions for each detection. If there is a data gap is longer than 1 year, we add a Heaviside step function located in the center of the gap. Figure 4.5 shows the remodeled transient velocity field, which is defined as the sum of all the step functions divided by the total duration of the time series; this represents the contribution of all the transient events to the overall time series trend. Overall, our result looks very similar to the RSI result, except some larger velocity toward to the trench at the stations around latitude of 48o and 40o. This is likely because Crowell et al. (2016) 70 had used a larger detection threshold of detection based on just white noise, so that some smaller sizes of transients may were missed in their study. Figure 4.5. Transient velocity field, remodeled based on transient detections with probability greater than 0.5, in the Cascadia from 2005 to 2016 from the results of (left) our ML model, (middle) the RSI model and (right) their difference. 4.4.3. Testing on additional data that are “unseen” in the training data In this section, we test how the model performs on some additional data under conditions that are different from the training data set. These tests could help us understand what pattern has been essentially learnt in the ML model. First, we apply the model to the site AB48, which is located in Port Alexander, east of the Gulf of Alaska. This station contains the coseismic and postseismic deformation from the 2013 Mw7.5 Craig earthquake. The north component of the time series and the model prediction are 71 shown in Figure 4.6a. The coseismic offset is identified as a very short duration transient. The earthquake, as a special type of transient deformation, is super short and the offset size is greater than 50 mm. Both are outside the range of the transient parameters used to generate the training samples, but we can see the predictions identifies the earthquake and its correct timing. The postseismic deformation that follows the earthquake is not detected as a transient, likely because its timescale of variation is much longer than that of the training data. Second, we test the site AC03 as an example that contains much longer transients. AC03 is located in Cook Inlet close to Homer, Alaska, and includes deformation due to multi-year slow slip events (Li et al., 2016). Figure 4.6b shows the east component of the time series and the model prediction for AC03. The model is not able to do the correct prediction using the original time series. However, the model does the correct prediction after we uniformly decimate the time series by selecting one out of every 20 epochs and then use the decimated time series as the input for prediction; the resampling effectively compresses the time series in the time domain. From the perspective of the model, the duration of the transient reduces from about 600 days to 30 days which lays inside the training range. Because our model is trained without using any input information about time, extending or compressing the time series in the time domain can allow the model to detect transients of various durations. Because our model is trained without using any input information about time, it has learned the pattern of transient motion, and this also suggests that long-duration transients are similar enough to short-duration transients except for timescale. Site ATW2 is located in Upper Cook Inlet, Alaska, and contains even longer transient of about 5 years from 2009 to 2014. Again, we resample the time series by selecting one of every other 60th epochs for the detection. The model successfully pick out the times of apparently accelerated deformation during this long transient (Figure 4.6c). Thus, resampling data to different degrees of 72 time compression may be an effective approach for detecting transients of varying lengths. Alternatively, this result suggests that it should be possible to train a model for a wider range of transient event durations than used here, although that would require much more computational time as the input training series would also have to be much longer. 4.5. Conclusion We present a simple but efficient ML method for transient detection based on the LSTM layer, a type of RNN. Our model takes each data point in a geodetic time series as input and directly output its transient probability in a sequential way. Through synthetic tests, we have shown that that our ML model is able to detect transients of smaller sizes compared with the RSI method under the same noise level, benefiting from the feature extraction with the ML layers. We later apply our model to a real-world dataset for detecting the slow slip events in the Cascadia region. Our detection result is overall consistent with the previous studies, which validates the model performance in a real study case as evidenced in the synthetic tests. The model successfully detects coseismic offsets but does not detect postseismic transients. By resampling/decimating a time series to a lower data rate to effectively compress time, the model is able to recognize long-duration transients even when it has been trained only on short-duration transients. 73 Figure 4.6. Test results for AB48, AC03 and ATW2. For AC03, the gray line in (b) is the prediction results using the original time series as input and the yellow line shows the results using the resampled time series. 74 Chapter 5: Conclusion To fully utilize InSAR time series with complementary GPS time series, I have introduced an improved inversion filter for modeling the volcano deformation that is able to combine GPS and InSAR. Using this improved inversion filter, I model the post-eruptive deformation at Okmok volcano, Alaska, since 2008. By combining GPS and InSAR timeseries observations, I determine the structure of the magma plumbing system should be a combination of a Mogi source and a shallow sill to account for the volcanic deformation. While the Mogi source is close to the general pre-eruptive Mogi source that should represent the long-lived magma body that already existed since before the 1997 eruption, the shallow sill is likely an intrusion into a shallow space that was opened when an older, more evolved magma body was partially drained during the 2008 eruption. This work demonstrates that combining GPS and InSAR is crucial to understand the shallow magma plumbing system at Okmok because the GPS data and the InSAR data individually suggest two distinct sources and the final “Mogi+sill” model cannot be determined without using both data sets. The inversion filter is also a very efficient tool to recover the volume change history. In total, five inflation episodes are identified during the post-eruptive period at Okmok, each with the form of an exponential decay. The total volume increase from the end of the 2008 eruption to 2019 is at least 50% of the 2008 co-eruptive volume decrease. Chapter 3 presents another work using the inversion filter. The location and the time- dependent volume change of the shallow magma chamber are estimated for all three volcanoes in east central Alaska Aleutian arc. The ~25 years volume change history shows two exponentially decaying pulses followed by a nearly linear inflation period at Westdahl, multiple episodic inflation events at Akutan, and transient volcanic inflation and long-term deflation signals at Makushin. These results indicate that the magma supply dynamics of the three volcanoes are 75 diverse even though they are close in distance. In addition, an accurate secular velocity is separated from the large volcanic deformation for the study region using the inversion filter. A small (~2 mm/year) northward component is found for the volcanoes relative to the overall block velocities estimated by the block models. This small signal is inferred to be due to the elastic strain from slip deficit on the megathrust, and two coupling models are tested: 1) a shallow (0–15 km) nearly fully locked zone, 2) a low slip deficit rate (~12% of the plate motion) on only a deeper locked zone at 35–50 km depth. Both models fit the data equally well, but the former model that was favored by the tsunami simulation results (Nicolsky et al., 2016). A single exponential inflation can be explained by a simple hydraulic model made up by a single shallow magmatic source fed by an infinite deep magma reservoir (which means constant influx) (e.g. Pinel et al., 2010), with an initial pressure difference between the two. The time- dependent volume change of all four volcanoes studied here show multiple inflation events. These results cannot be explained by a single pressure pulse or a steady flux from below, but require multiple pulses of magma from depth. More specifically, the pressure of the deep source should increase in episodic pulses so that each new inflation episode occurs as a result of episodic magma supply into this deep source. We also see that the timescale of the inflation events at Okmok decreases over time. This could be explained as the influx from deep to shallow sources increase with time. A more thorough model is worthy of future work to fully explain the continuously changing magma supplies. The volume change deficit is considered to be an indicator for the timing of the next eruption. The capacity of the magma storage system is assumed within a certain range, and an eruption will likely occur if that capacity has been reached since the last eruption. At Okmok, in the previous eruptive cycle, the inter-eruptive volume increase from 1997 to 2008 was 85-100% 76 of the 1997 eruption volume (Lu et al., 2010). During the current cycle, an estimation of 0.068 km3 of magma has been added to the shallow magma storage from 2008 to 2019, which is 50-60% of the volume withdrawn during the 2008 eruption. The hypothesis of volume deficit can be further validated using the inversion filter presented here and future observations. In Chapter 4, I present a simple but efficient ML model for transient detection based on a recurrent neural network with LSTM cells. Unlike previous studies, this new model sequentially takes each data point in a time series as input only once, and directly output the transient probability instead of adopting a sliding time window strategy. The model is trained on synthetic data of fixed short length but can be applied to timeseries of any lengths for prediction. Synthetic tests have shown that our ML model is able to detect transients of smaller sizes compared with the RSI method under the same noise level. I then apply the model to a real-world dataset for detecting the slow slip events in the Cascadia region. The detection result is overall consistent with the previous studies. This validates the model performance in a real study case as evidenced in the synthetic tests. The additional tests against completely unseen time series shows that the ML model does generalize a pattern that is able to detect transients of varying characteristics. The work of Chapter 4 can be directly expanded by training using real data instead of synthetic data, although it may be difficult to obtain a real data set that is well labeled in a standard way. Chapter 4 gives an application of the proposed ML method in geodetic data, but the method can be easily transferred to studies on other types of time series observations with few modifications. It would be interesting to see how the proposed method performs in the application of, for example, seismic phase picking as the seismic data is considered to be a lager in size and have higher frequencies compared with the geodetic data. 77 BIBLIOGRAPHY 78 BIBLIOGRAPHY Albright, J. A., Gregg, P. M., Lu, Z., & Freymueller, J. T. (2019). Hindcasting Magma Reservoir Stability Preceding the 2008 Eruption of Okmok, Alaska. Geophysical Research Letters, 46(15), 8801–8808. https://doi.org/https://doi.org/10.1029/2019GL083395 Anantrasirichai, N., Biggs, J., Albino, F., Hill, P., & Bull, D. (2018). Application of Machine Learning to Classification of Volcanic Deformation in Routinely Generated InSAR Data. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2018JB015911 Anantrasirichai, N., Biggs, J., Albino, F., & Bull, D. (2019). A deep learning approach to detecting volcano deformation from satellite imagery using synthetic datasets. Remote Sensing of Environment. https://doi.org/10.1016/j.rse.2019.04.032 Anderson, K., & Segall, P. (2011). Physics-based models of ground deformation and extrusion rate at effusively erupting volcanoes. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2010JB007939 Avé Lallemant, H. G., & Oldow, J. S. (2000). Active displacement partitioning and arc-parallel extension of the Aleutian volcanic arc based on Global Positioning System geodesy and kinematic analysis. Geology. https://doi.org/10.1130/0091- 7613(2000)28<739:ADPAAE>2.0.CO;2 Bartlow, N. M., Wallace, L. M., Beavan, R. J., Bannister, S., & Segall, P. (2014). Time- dependent modeling of slow slip events and associated seismicity and tremor at the Hikurangi subduction zone, New Zealand. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1002/2013JB010609 Bekaert, D. P. S., Segall, P., Wright, T. J., & Hooper, A. J. (2016). A Network Inversion Filter combining GNSS and InSAR for tectonic slip modeling. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1002/2015JB012638 Biggs, J., Lu, Z., Fournier, T., & Freymueller, J. T. (2010). Magma flux at Okmok Volcano, Alaska, from a joint inversion of continuous GPS, campaign GPS, and interferometric synthetic aperture radar. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2010JB007577 Blewitt, G., Hammond, W., & Kreemer, C. (2018). Harnessing the GPS Data Explosion for Interdisciplinary Science. Eos, 99. https://doi.org/10.1029/2018eo104623 Bradley, A. P. (1997). The use of the area under the ROC curve in the evaluation of machine learning algorithms. Pattern Recognition, 30(7), 1145–1159. https://doi.org/https://doi.org/10.1016/S0031-3203(96)00142-2 Bürgmann, R., & Chadwell, D. (2014). Seafloor Geodesy. Annual Review of Earth and Planetary Sciences, 42(1), 509–534. https://doi.org/10.1146/annurev-earth-060313-054953 79 Cervelli, P. F. (2013). Analytical Expressions for Deformation from an Arbitrarily Oriented Spheroid in a Half-Space. Eos, Transactions American Geophysical Union. Cervelli, Peter F., & Miklius, A. (2003). The shallow magmatic system of Kīlauea Volcano. US Geological Survey Professional Paper. Cross, R. S., & Freymueller, J. T. (2008). Evidence for and implications of a Bering plate based on geodetic measurements from the Aleutians and western Alaska. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2007JB005136 Crowell, B. W., Bock, Y., & Liu, Z. (2016). Single-station automated detection of transient deformation in GPS time series with the relative strength index: A case study of Cascadian slow slip. Journal of Geophysical Research: Solid Earth, 121(12), 9077–9094. https://doi.org/10.1002/2016JB013542 DeGrandpre, K., Wang, T., Lu, Z., & Freymueller, J. T. (2017). Episodic inflation and complex surface deformation of Akutan volcano, Alaska revealed from GPS time-series. Journal of Volcanology and Geothermal Research. https://doi.org/10.1016/j.jvolgeores.2017.10.003 Delaney, P. T., & McTigue, D. F. (1994). Volume of magma accumulation or withdrawal estimated from surface uplift or subsidence, with application to the 1960 collapse of Kilauea volcano. Bulletin of Volcanology, 56(6–7), 417–424. https://doi.org/10.1007/BF00302823 Delouis, B., Giardini, D., Lundgren, P., & Salichon, J. (2002). Joint inversion of InSAR, GPS, teleseismic, and strong-motion data for the spatial and temporal distribution of earthquake slip: Application to the 1999 İzmit mainshock. Bulletin of the Seismological Society of America. https://doi.org/10.1785/0120000806 Dong, D., Fang, P., Bock, Y., Webb, F., Prawirodirdjo, L., Kedar, S., & Jamason, P. (2006). Spatiotemporal filtering using principal component analysis and Karhunen-Loeve expansion approaches for regional GPS network analysis. Journal of Geophysical Research: Solid Earth, 111(B3). https://doi.org/https://doi.org/10.1029/2005JB003806 Dzurisin, D., & Poland, M. P. (2019). Magma supply to Kīlauea Volcano, Hawai‘i, from inception to now: Historical perspective, current state of knowledge, and future challenges. In Field Volcanology: A Tribute to the Distinguished Career of Don Swanson. Geological Society of America. https://doi.org/10.1130/2018.2538(12) Elliott, J., & Freymueller, J. T. (2020). A Block Model of Present‐Day Kinematics of Alaska and Western Canada. Journal of Geophysical Research: Solid Earth, n/a(n/a), e53951. https://doi.org/10.1029/2019JB018378 Evensen, G. (1994). Sequential data assimilation with a nonlinear quasi-geostrophic model using Monte Carlo methods to forecast error statistics. Journal of Geophysical Research: Oceans, 99(C5), 10143–10162. https://doi.org/https://doi.org/10.1029/94JC00572 Fialko, Y., Khazan, Y., & Simons, M. (2001). Deformation due to a pressurized horizontal circular crack in an elastic half-space, with applications to volcano geodesy. Geophysical 80 Journal International. https://doi.org/10.1046/j.1365-246X.2001.00452.x Fournier, T., Freymueller, J., & Cervelli, P. (2009). Tracking magma volume recovery at Okmok volcano using GPS and an unscented Kalman filter. Journal of Geophysical Research, 114(B2), B02405. https://doi.org/10.1029/2008JB005837 Freymueller, J. T., & Kaufman, A. M. (2010). Changes in the magma system during the 2008 eruption of Okmok volcano, Alaska, based on GPS measurements. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2010JB007716 Fu, Y., & Freymueller, J. T. (2012). Seasonal and long-term vertical deformation in the Nepal Himalaya constrained by GPS and GRACE measurements. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2011JB008925 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. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1002/2014JB011384 Gregg, P. M., & Pettijohn, J. C. (2016). A multi-data stream assimilation framework for the assessment of volcanic unrest. Journal of Volcanology and Geothermal Research. https://doi.org/10.1016/j.jvolgeores.2015.11.008 He, B., Wei, M., Watts, D. R., & Shen, Y. (2020). Detecting Slow Slip Events From Seafloor Pressure Data Using Machine Learning. Geophysical Research Letters, 47(11). https://doi.org/10.1029/2020GL087579 Henderson, S. T., & Pritchard, M. E. (2013). Decadal volcanic deformation in the Central Andes Volcanic Zone revealed by InSAR time series. Geochemistry, Geophysics, Geosystems, 14(5), 1358–1374. https://doi.org/10.1002/ggge.20074 Hochreiter, S., & Schmidhuber, J. (1997). Long Short-Term Memory. Neural Computation, 9(8), 1735–1780. https://doi.org/10.1162/neco.1997.9.8.1735 Ji, K. H., & Herring, T. A. (2011). Transient signal detection using GPS measurements: Transient inflation at Akutan volcano, Alaska, during early 2008. Geophysical Research Letters. https://doi.org/10.1029/2011GL046904 Ji, K. H., Yun, S. H., & Rim, H. (2017). Episodic inflation events at Akutan Volcano, Alaska, during 2005–2017. Geophysical Research Letters. https://doi.org/10.1002/2017GL074626 Julier, S. J., & Uhlmann, J. K. (2004). Unscented filtering and nonlinear estimation. In Proceedings of the IEEE. https://doi.org/10.1109/JPROC.2003.823141 Kingma, D. P., & Ba, J. L. (2015). Adam: A method for stochastic optimization. 3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings. 81 Koulakov, I., Abkadyrov, I., Al Arifi, N., Deev, E., Droznina, S., Gordeev, E. I., et al. (2017). Three different types of plumbing system beneath the neighboring active volcanoes of Tolbachik, Bezymianny, and Klyuchevskoy in Kamchatka. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1002/2017JB014082 Langbein, J. (2008). Noise in GPS displacement measurements from Southern California and Southern Nevada. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2007JB005247 Larsen, J. F., Neal, C. A., Schaefer, J. R., Kaufman, A. M., & Lu, Z. (2015). The 2008 phreatomagmatic eruption of Okmok Volcano, Aleutian Islands, Alaska: Chronology, deposits, and landform changes. https://doi.org/10.14509/29405 Larsen, Jessica F., Śliwiński, M. G., Nye, C., Cameron, C., & Schaefer, J. R. (2013). The 2008 eruption of Okmok Volcano, Alaska: Petrological and geochemical constraints on the subsurface magma plumbing system. Journal of Volcanology and Geothermal Research. https://doi.org/10.1016/j.jvolgeores.2013.07.003 Lengliné, O., Marsan, D., Got, J. L., Pinel, V., Ferrazzini, V., & Okubo, P. G. (2008). Seismicity and deformation induced by magma accumulation at three basaltic volcanoes. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2008JB005937 Li, S., & Freymueller, J. T. (2018). Spatial Variation of Slip Behavior Beneath the Alaska Peninsula Along Alaska-Aleutian Subduction Zone. Geophysical Research Letters. https://doi.org/10.1002/2017GL076761 Li, S., Freymueller, J., & McCaffrey, R. (2016). Slow slip events and time‐dependent variations in locking beneath Lower Cook Inlet of the Alaska‐Aleutian subduction zone. Journal of Geophysical Research: Solid Earth, 121(2), 1060–1079. https://doi.org/10.1002/2015JB012491 Lohman, R. B., & Simons, M. (2005). Some thoughts on the use of InSAR data to constrain models of surface deformation: Noise structure and data downsampling. Geochemistry, Geophysics, Geosystems. https://doi.org/10.1029/2004GC000841 Lu, Z., & Dzurisin, D. (2010). Ground surface deformation patterns, magma supply, and magma storage at Okmok volcano, Alaska, from InSAR analysis: 2. Coeruptive deflation, July– August 2008. Journal of Geophysical Research: Solid Earth, 115(B5). https://doi.org/https://doi.org/10.1029/2009JB006970 Lu, Z., & Dzurisin, D. (2014). InSAR Imaging of Aleutian Volcanoes. InSAR Imaging of Aleutian Volcanoes. Berlin, Heidelberg: Springer Berlin Heidelberg. https://doi.org/10.1007/978-3- 642-00348-6 Lu, Z., Power, J. A., McConnell, V. S., Wicks, C., & Dzurisin, D. (2002). Preeruptive inflation and surface interferometric coherence characteristics revealed by satellite radar interferometry at Makushin Volcano, Alaska: 1993-2000. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2001jb000970 82 Lu, Z., Masterlark, T., Dzurisin, D., Rykhus, R., & Wicks, C. (2003). Magma supply dynamics at Westdahl volcano, Alaska, modeled from satellite radar interferometry. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2002jb002311 Lu, Z., Masterlark, T., & Dzurisin, D. (2005). Interferometric synthetic aperture radar study of Okmok volcano, Alaska, 1992-2003: Magma supply dynamics and postemplacement lava flow deformation. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2004JB003148 Lu, Z., Dzurisin, D., Wicks, C., Power, J., Kwoun, O., & Rykhus, R. (2007). Diverse deformation patterns of Aleutian Volcanoes from satellite Interferometric Synthetic Aperture Radar (InSAR). In Volcanism and Subduction: The Kamchatka Region (pp. 249– 261). https://doi.org/10.1029/172GM18 Lu, Z., Dzurisin, D., Biggs, J., Wicks, C., & McNutt, S. (2010). Ground surface deformation patterns, magma supply, and magma storage at Okmok volcano, Alaska, from InSAR analysis: 1. Intereruption deformation, 19970-2008. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2009JB006969 Mann, D., & Freymueller, J. (2003). Volcanic and tectonic deformation on Unimak Island in the Aleutian Arc, Alaska. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2002jb001925 Mann, D., Freymueller, J., & Lu, Z. (2002). Deformation associated with the 1997 eruption of Okmok volcano, Alaska. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/2001jb000163 McGuire, J. J., & Segall, P. (2003). Imaging of aseismic fault slip transients recorded by dense geodetic networks. Geophysical Journal International. https://doi.org/10.1111/j.1365- 246X.2003.02022.x Melbourne, T. I., Szeliga, W. M., Miller, M. M., & Santillan, V. M. (2005). Extent and duration of the 2003 Cascadia slow earthquake. Geophysical Research Letters, 32(4). https://doi.org/https://doi.org/10.1029/2004GL021790 Miller, T. P., McGimsey, R. G., Richter, D. H., Riehle, J. R., Nye, C. J., Yount, M. E., & Dumoulin, J. A. (1998). Catalog of the historically active volcanoes of Alaska. Open-File Report. https://doi.org/10.3133/ofr98582 Mogi, K. (1958). Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Bulletin of the Earthquake Research Institute. https://doi.org/10.1016/j.epsl.2004.04.016 Nicolsky, D. J., Freymueller, J. T., Witter, R. C., Suleimani, E. N., & Koehler, R. D. (2016). Evidence for shallow megathrust slip across the Unalaska seismic gap during the great 1957 Andreanof Islands earthquake, eastern Aleutian Islands, Alaska. Geophysical Research Letters, 43(19), 10,328-10,337. https://doi.org/10.1002/2016GL070704 83 Pinel, V., Jaupart, C., & Albino, F. (2010). On the relationship between cycles of eruptive activity and growth of a volcanic edifice. Journal of Volcanology and Geothermal Research. https://doi.org/10.1016/j.jvolgeores.2010.05.006 Qu, F., Lu, Z., Poland, M., Freymueller, J., Zhang, Q., & Jung, H.-S. (2015). Post-Eruptive Inflation of Okmok Volcano, Alaska, from InSAR, 2008–2014. Remote Sensing. https://doi.org/10.3390/rs71215839 Riel, B., Simons, M., Agram, P., & Zhan, Z. (2014). Detecting transient signals in geodetic time series using sparse estimation techniques. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1002/2014JB011077 Ross, Z. E., Meier, M.-A., & Hauksson, E. (2018). P Wave Arrival Picking and First-Motion Polarity Determination With Deep Learning. Journal of Geophysical Research: Solid Earth, 123(6), 5120–5129. https://doi.org/10.1029/2017JB015251 Schmidt, D. A., & Gao, H. (2010). Source parameters and time-dependent slip distributions of slow slip events on the Cascadia subduction zone from 1998 to 2008. Journal of Geophysical Research, 115, B00A18. https://doi.org/10.1029/2008JB006045 Segall, P., & Matthews, M. (1997). Time dependent inversion of geodetic data. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1029/97jb01795 Sigmundsson, F., Hreinsdóttir, S., Hooper, A., Árnadóttir, T., Pedersen, R., Roberts, M. J., et al. (2010). Intrusion triggering of the 2010 Eyjafjallajökull explosive eruption. Nature, 468(7322), 426–430. https://doi.org/10.1038/nature09558 Simons, M., Fialko, Y., & Rivera, L. (2002). Coseismic Deformation from the 1999 Mw 7.1 Hector Mine, California, Earthquake as Inferred from InSAR and GPS Observations. Bulletin of the Seismological Society of America, 92(4), 1390–1402. https://doi.org/10.1785/0120000933 Walwer, D., Ghil, M., & Calais, E. (2019). Oscillatory nature of the Okmok volcano’s deformation. Earth and Planetary Science Letters. https://doi.org/10.1016/j.epsl.2018.10.033 Walwer, Damian, Calais, E., & Ghil, M. (2016). Data-adaptive detection of transient deformation in geodetic networks. Journal of Geophysical Research: Solid Earth. https://doi.org/10.1002/2015JB012424 Wan, E. A., & Van Der Merwe, R. (2000). The unscented Kalman filter for nonlinear estimation. In IEEE 2000 Adaptive Systems for Signal Processing, Communications, and Control Symposium, AS-SPCC 2000. https://doi.org/10.1109/ASSPCC.2000.882463 Wang, T., DeGrandpre, K., Lu, Z., & Freymueller, J. T. (2018). Complex surface deformation of Akutan volcano, Alaska revealed from InSAR time series. International Journal of Applied Earth Observation and Geoinformation. https://doi.org/10.1016/j.jag.2017.09.001 84 Xu, C., Ding, K., Cai, J., & Grafarend, E. W. (2009). Methods of determining weight scaling factors for geodetic-geophysical joint inversion. Journal of Geodynamics. https://doi.org/10.1016/j.jog.2008.06.005 Xu, C., Liu, Y., Wen, Y., & Wang, R. (2010). Coseismic slip distribution of the 2008 mw 7.9 Wenchuan earthquake from joint inversion of GPS and InSAR data. Bulletin of the Seismological Society of America. https://doi.org/10.1785/0120090253 Xue, X., Freymueller, J., & Lu, Z. (2020). Modeling the Posteruptive Deformation at Okmok Based on the GPS and InSAR Time Series: Changes in the Shallow Magma Storage System. Journal of Geophysical Research: Solid Earth, 125(2), e2019JB017801. https://doi.org/10.1029/2019JB017801 Yang, X. M., Davis, P. M., & Dieterich, J. H. (1988). Deformation from inflation of a dipping finite prolate spheroid in an elastic half-space as a model for volcanic stressing. Journal of Geophysical Research. https://doi.org/10.1029/JB093iB05p04249 Zhan, Y., & Gregg, P. M. (2017). Data assimilation strategies for volcano geodesy. Journal of Volcanology and Geothermal Research. https://doi.org/10.1016/j.jvolgeores.2017.02.015 Zhu, W., & Beroza, G. C. (2018). PhaseNet: A Deep-Neural-Network-Based Seismic Arrival Time Picking Method. Geophysical Journal International. https://doi.org/10.1093/gji/ggy423 85