Improved Real-Time Natural Hazard Monitoring Using Automated DInSAR Time Series

: As part of the collaborative GeoSciFramework project, we are establising a monitoring system for the Yellowstone volcanic area that integrates multiple geodetic and seismic data sets into an advanced cyber-infrastructure framework that will enable real-time streaming data analytics and machine learning and allow us to better characterize associated long- and short-term hazards. The goal is to continuously ingest both remote sensing (GNSS, DInSAR) and ground-based (seismic, thermal and gas observations, strainmeter, tiltmeter and gravity measurements) data and query and analyse them in near-real time. In this study, we focus on DInSAR data processing and the effects from using various atmospheric corrections and real-time orbits on the automated processing and results. We ﬁnd that the atmospheric correction provided by the European Centre for Medium-Range Weather Forecasts (ECMWF) is currently the most optimal for automated DInSAR processing and that the use of real-time orbits is sufﬁcient for the early-warning application in question. We show analysis of atmospheric corrections and using real-time orbits in a test case over the Kilauea volcanic area in Hawaii. Finally, using these ﬁndings, we present results of displacement time series in the Yellowstone area between May 2018 and October 2019, which are in good agreement with GNSS data where available. These results will contribute to a baseline model that will be the basis of a future early-warning system that will be continuously updated with new DInSAR data acquisitions.


Introduction
Natural catastrophes occur at a variety of spatial and temporal scales. In particular, solid earth hazards, such as large earthquakes and volcanic eruptions, often have very long inter-event times and this makes it difficult to forecast their behavior. Today, we monitor a significant number of active volcanoes around the world, collecting geodetic data in the form of GNSS and SAR data, in conjunction with ground-based measurements from strainmeters, tiltmeters and gravity instruments with the goal of providing early warning and reducing risk and losses. In addition, in many areas, we also collect seismic data, thermal measurements and/or perform sampling of local gas emissions. While the term "big data" is widely used today, it is not simply about volume, but also about the variability in data types, the latency with which data is generated and transmitted, and the veracity, or accuracy of the models that ingest data [1]. Natural volcano observatories that collect large volumes of variable data over long time periods provide a unique opportunity to develop a framework that integrates modern data sources and advanced computational algorithms. This is essential to both early detection of volcanic activity and forecasting eruptions on intermediate-and short-term time scales.
Modern forecasting tools for volcanic eruptions primarily rely on some combination of probabilistic studies of geologic records of activity and eruption and comprehensive monitoring networks that incorporate seismic, geodetic, thermal and gas measurements into modern computational tools [2][3][4][5][6][7][8][9]. In particular, recent developments in differential interferometric synthetic aperture radar (DInSAR) and Global Navigation Satellite System (GNSS) geodes have promoted significant advances in monitoring of volcanic hazards, providing high-resolution ground deformation data with dense spatiotemporal coverage. DInSAR can be used to map ground deformation with high spatial resolution and sub-centimeter precision over large areas, and is a suitable tool for ground deformation monitoring of active volcanic regions [10][11][12][13][14][15]. The unprecedented volumes of SAR data provided by modern satellites such as Sentinel-1A/B from the European Space Agency (ESA), ALOS-2 from the Japan Aerospace Exploration Agency (JAXA) and the upcoming NISAR (NASA-ISRO Synthetic Aperture Radar) mission from the US National Aeronautics and Space Administration (NASA) and the Indian Space Research Organization present both an important opportunity and a significant challenge for observation of the large number of active volcanoes around the world that are not currently monitored by volcano observatories [6,8] A number of research groups around the world are developing high-performance computing and cloud-based methods for rapid processing of both SAR and DInSAR images. These include, among others, the Advanced Rapid Imaging and Analysis (ARIA) project [16][17][18], the Alaskan Satellite Facility (ASF) [19][20][21], the Centre for the Observation and Modelling of Earthquakes, Volcanoes and Tectonics (COMET) [22,23] and the UNAVCO GeoSciFramework project [24]. The natural extension of these efforts focuses on the longerterm monitoring of specific regions, such as the national ground motion service for Norway (Annual Report, Geological Survey of Norway, https://www.ngu.no/en/topic/insarnorway accessed on 23 February 2021), and specific phenomena, including volcanic activity in South America [25].
Acquisition and processing of satellite SAR data and GNSS data (sampling over time frames of days-to-years) can be used to study strain accumulation in Earth's crust, allowing us to model and invert for the underlying physical cause of that strain, including changes in magma composition, storage and movement. Recent developments in remote sensing make it feasible to continuously monitor and predict natural hazards. In particular, with the Sentinel missions of the European Space Agency (ESA), remote sensing data is becoming freely available to researchers everywhere [26]. The Sentinel-1 missions are providing SAR measurements for the entire globe with an unprecedented repeatability of 12 days. When combining Sentinel 1A and Sentinel 1B acquisitions, the repeatability becomes 6 days. These repeatedly acquired SAR data from a single sensor can be used to obtain the Line-of-Sight (LOS) time series of ground deformation using DInSAR methods [14,[27][28][29][30].
As part of the GeoSciFramework project [31] (https://www.unavco.org/projects/ major-projects/earthcube/geosciframework/geosciframework.html accessed on 23 February 2021), we are focusing on modeling and monitoring ground motion in the Yellowstone volcanic area. The ultimate goal of the project is to continuously monitor this area using newly acquired DInSAR and GNSS data and, in conjunction with thermal and gas observations and strainmeter, tiltmeter and gravity measurements, integrate that data into machine learning approaches to identify precursory patterns and inform eruption forecasts. Here, we process DInSAR images for Yellowstone into time series, then include them in the time series stream integrated with GNSS data; this provides 3D surface motions, improves the accuracy of the SAR-derived time series and derives high-resolution time series maps, complete with the corresponding horizontal and vertical errors [32].
In this paper, we discuss details of the DInSAR data processing for this project with specific focus on the effect of using real time orbits and standard atmospheric corrections in the automated processing. When monitoring a region of potential natural hazard such as a volcano, it is necessary to use real-time orbits when processing DInSAR data, as these orbits are available weeks earlier than the precise orbits. Therefore, it is important to understand the potential limitations and accuracy of DInSAR data processed with these orbits. DInSAR data can be noisy due to atmospheric effects and it is essential that we apply atmospheric corrections during the processing. In this paper, we explore the use of atmospheric corrections from two providers: European Centre for Medium-Range Weather Forecasts (ECMWF), a correction that is widely used [33], and the Generic Atmospheric Correction Online Service (GACOS), which integrates ECMWF and continuous GNSS tropospheric delay estimates [34]. While GACOS might become the state-of-the-art atmospheric correction in the coming years, it is not yet implemented for automated processing, and for the near-real time application that we discuss in this paper, it could be more beneficial to use the ECMWF correction that can be applied automatically during the processing flow.
Because the Yellowstone region is challenging for DInSAR data processing due to its dense vegetation and regular snow coverage, we tested the accuracy of DInSAR data with different orbits over the Big Island of Hawaii. Hawaii is a simpler case for DInSAR data processing, is well-instrumented over the entire island, and there are other studies available for comparison with our results. We apply these results to an initial study of Yellowstone caldera [35].

The Yellowstone and Hawaiian Volcanic Areas
The Yellowstone caldera is a long-lived, silicic volcanic system located on the eastern edge of the Basin and Range province. Fed by a deep-seated mantle plume that is responsible for its long-term activity, there have been three caldera-forming eruptions, at 2.0, 1.3 and 0.64 Ma BP, at this intraplate hotspot. That most recent eruption released almost 1000 km 3 of material and collapsed into a 45 by 70 km caldera [36][37][38][39]. Considered one of the most potentially dangerous volcanoes in the world [40], the caldera is characterized by very high heat flow, widespread seismicity, extensive hydrothermal activity, and ongoing episodes of surface uplift and subsidence [38]. For example, since 2004, there have been two episodes of surface deformation with uplift rates of as much as 7 cm/year. Modeling of that deformation identified an expanding sill at 7-10 km in depth [41].
Spatial and temporal variations of Yellowstone ground movement are correlated with changes in seismic and hydrothermal activity in and around the caldera [38,[41][42][43]. Because of the dynamic nature of the caldera and the significant volcanic and seismic hazard, there is a substantial monitoring effort at Yellowstone, including continuous GNSS, borehole strain, and seismic networks. Figure 1 shows the Sentinel-1 frames that cover the Yellowstone area.
The island of Hawaii is formed by five volcanoes, among them the youngest and most active is the Kilauea volcano. An active shield volcano that was created as the Pacific plate moved over the Hawaiian hotspot [44], Kilauea was formed at least 300,000 years ago, reaching above sea-level about 100,000 years ago [45]. The volcano has been continuously active in modern history, producing a stream of explosive and effusive eruptions [46].
The island and Kilauea volcano has been extensively studied and is continuously monitored by the USGS Hawaiian Volcano Observatory. The monitoring networks consist of more than 100 field stations with instruments that record and measure earthquakes, ground movement, volcanic gases, sound waves, lava advancement, inferred magma volume below ground, and visual changes in eruptive activity. In addition, remote sensing, especially DInSAR data, can provide an additional large spatial and temporal-scale overview [47,48].
In Figure 2, we show the descending frame of the Sentinel-1 DInSAR data coverage from the data we use in this paper. This is not a complete image of the DInSAR coverage over the Hawaii area.

DInSAR Data
We acquired Sentinel-1 data through the Alaskan Satellite Facility (ASF) and processed them using the GMT5SAR DInSAR data processing package [49]. We generated time series of displacement using the GIAnT toolbox [50]. We processed 43 scenes of Sentinel-1 DInSAR data between 4 January 2017 and 16 June 2018 of the descending frame of Path 87, Frame 527 ( Figure 2). We generated every possible combination of interferograms with less than 50 days of temporal baseline. Because the spatial repeatability of the Sentinel-1 acquisitions does not vary significantly over the 18 months of the period that we were interested in, we did not impose any spatial baseline limits.
We processed all available Sentinel-1 SAR data covering the Yellowstone volcanic area. Four frames are fully or partially covering the Yellowstone National Park and the surrounding areas ( Figure 1). We used the GMT5SAR [49] toolbox for DInSAR processing and the GIAnT [50] toolbox for displacement time series calculation. These were carried out using the Summit high-performance research computing facilities at the University of Colorado, Boulder.
We used all available SAR acquisitions between January 2017 and December 2018, and formed all possible interferograms with a maximum temporal baseline of 50 days. We did not impose any spatial baseline limit, as the Sentinel-1 mission is designed to have a great repeatability in terms of spatial coverage of the same frame. Details of data processing of each frame are shown in Table 1.
Details of data processed in this study for both the Hawaii and Yellowstone areas are shown in Table 1.

Atmospheric Correction and Displacement Time Series
The atmosphere varies greatly both in space and time, and radar signals that propagate through the various layers of the atmosphere are significantly affected. The atmosphere can introduce errors of over ten centimeters in the measured ground deformation [51]. Atmospheric effects can be mitigated based on external data, such as meteorological observations [52,53], GNSS observations [54,55], high-resolution meteorological models [56], or water vapor (MODIS [57] or MERIS [55]) data. Correction can also be applied based on techniques developed for DInSAR data analysis, such as correlation analysis [58,59], pair-wise logic [60], PSInSAR technique [61,62] or stacking [61,63].
Various atmospheric corrections are available. In this paper, we explore the application of two of them, both based on external datasets, such as meteorological observations and GNSS data.
We assess, using the ECMWF [33] and the GACOS [34], atmospheric corrections based on their effectiveness and availability. The ECMWF (European Centre for Medium-Range Weather Forecasts) is based on an 80 km resolution weather model, and the GACOS correction is incorporating the ECMWF weather data with continuous GNSS measurements.
To calculate displacement time series, we use the Generic InSAR Analysis Toolbox (GIAnT) software package [50]. GIAnT is used for rapid generation of time-series of surface displacement using DInSAR data based on either the Small BAseline Subset (SBAS) or the New Small BAseline Subset (NSBAS) methods.
The SBAS technique is appropriate to apply when the differential interferograms have a small spatial baseline [27], that is typically the case with the Sentinel-1 mission data. The observation equation of the SBAS method is Equation (1).
where Φ ij is the phase of the interferogram combining acquisitions i and j and δϕ n is the phase increment between acquisitions n and n + 1. The cumulative phase can be obtained by solving a linear least squares inversion problem, but it solves only for coherent pixels, and this can result in a small number of pixels in the final solution.
The NSBAS method uses the same principle as the SBAS method but it (1) estimates DEM errors as well and (2) accommodates missing observations, and therefore, has been optimized for the monitoring of transient ground motion that has small amplitude and takes place over large areas and in natural settings [64]. The observation equation becomes as in Equation (2): where f (∆t k ) is the deformation and eB k perp is the DEM error. We expect this method to be appropriate to employ in the long-term monitoring of the Yellowstone area, but plan to explore other methods before the system becomes operational in order to optimize the output in terms of accuracy and processing time. For example, a new technique, Multidimensional Small Baseline Subset (MSBAS), allows us to incorporate DInSAR results from different satellites and wavelengths into one time series [65].

Real-time and Precise Orbits
We use the same Kilauea dataset as described in Section 2.3. We applied de-ramping and ECMWF atmospheric corrections [33]. Using these interferograms, we computed time series of displacements between 4 January 2017 and 16 June 2018 with the NSBAS technique, which is appropriate for a selected combination of interferograms with small spatial baselines between the images [64].
To compare results with precise and real-time orbits, we carried out the described processing sequence twice. In the first case, we used precise orbits for all the scenes and in the second case, we replaced the precise orbits with real-time ones for the last five scenes (Figure 3). This simulates a scenario of processing DInSAR data in near-real time, when for the past 20 days, there are only real-time orbits available but for older scenes there are the calculated precise orbits available. When operational, we will ingest new data when acquired and use real-time orbits with the newest DInSAR scenes, while at the same time, we will reprocess older scenes for which precise orbits became available since the last processing. Figure 3. Acquisitions processed, blue and red dots mark which orbit were used when discussing precise orbit and realtime orbit processing, respectively. Figure 4 shows results of three cases of processed DInSAR data over the Hawaii area: no atmospheric correction applied, ECMWF atmospheric correction applied [33], and GACOS atmospheric correction applied [34]. We find that both atmospheric corrections help remove apparent uplift signal (red) at the highest elevations of the island and keep the subsidence signal (blue) close to the southeast shoreline that is related to the 2018 Kilauea eruptions. Applying the ECMWF correction yields a smoother image than the GACOS correction.  Figure 5 shows the locations of selected GNSS stations and the corresponding time series of relative displacements derived from the DInSAR dataset. We selected three GNSS stations on the island of Hawaii that offer a comprehensive analysis of the corrections applied. At PMAU station, there is no expected ground motion; however, the ECMWF and GACOS cumulative displacement results seem to differ. MOKP station, which is close to the summit of Mauna Loa, is far enough from the volcanic activity that only small ground motion can be expected. Thirdly, at PUOC station, where ground motion is expected due to volcanic activity in 2018. We selected a single pixel closest to the location of each of these GNSS stations. Since some smoothing is already applied in the NSBAS process, we did not average among a number of pixels around the GNSS location. The DInSAR time series start at zero, and the GNSS time series are shifted to be shown overlaying the DInSAR time series. These GNSS time series are acquired from the Nevada Geodetic Laboratory and they span from 2000 for MOKP and from 2006 for both PUOC and PMAU, up until the present day [66].

Atmospheric Correction
We find that both ECMWF and GACOS perform similarly in terms of fit to the GNSS data, and even when no correction is applied, the uncorrected DInSAR results approximately match the corrected ones. In addition, we calculated the fit between the presented DInSAR and GNSS time series using the Dynamic Time Warping method (DTW) [67,68]. This technique calculates the similarities of time series of different sampling and gives the Euclidean distance between the two time series [69], d, as shown in each panel in Figure 5. The smaller the distance, the more similar the time series are. While none of the tests is best in every case, at the location of the 2018 volcanic activities (PUOC GNSS station), the ECMWF correction gives the smallest d value. Future studies should apply similar tests to larger suites of data to better assess under what conditions each method is preferred.
For further use in the GeoSciFramework project, we chose to use the ECMWF correction, as it is possible to include this in automatic processing chains as opposed to GACOS, which requires an individual data request for each scene. The DTW comparison shows that both ECMWF and GACOS perform very similarly, and as GACOS is planned to be automatically available in the future, we will consider this comparison in greater depth before the Yellowstone monitoring system becomes operational.

Real-Time and Precise Orbits
After processing 43 DInSAR acquisitions, we compare the two sets of results: one where we only used precise orbits, and one in which we replaced the orbits of the last five acquisitions with real-time orbits. Individual interferograms show some pixel-to-pixel difference, but the broad image of ground deformation and its interpretation do not differ between the two cases ( Figure 6). These two interferograms show the deformation of the active area that corresponds to the rift eruption and collapse of Kilauea volcano during spring and summer of 2018 [70]. We find similarly non-significant differences when assessing the comparison between cumulative displacements as a result of the two scenarios ( Figure 7). Differences at time steps 1, 15, 30 and 42 are shown in Figure 8. The 42nd time step is the final one, also corresponding to the interferogram in Figure 6. Because the first 38 acquisitions are the same in both processing sequences, we do not expect any disparities in the first 37 timesteps. However, as we carried out the processing and time series generation separately, and as the solution to the displacement history is non-unique, we still see some pixel-topixel differences at earlier time steps. The difference between the two sets of processing become larger and larger as we assess later and later timesteps; however, we find that these differences are due to the non-uniqueness of the solution for the displacement time series rather than a consequence of the real-time orbits used. If we assess the last five time-steps, where in one case, we use precise and in the other case, real-time orbits, we do not observe a change in the rate of increasing differences.

Yellowstone Results
We acquired Sentinel-1 data between January 2017 and December 2018 for the frames covering the Yellowstone National park (Figure 1). In the case of frames 100-441 and 100-446, we merged and cropped away the northern half of 100-441 and the southern half of 100-446, and processed the mid-section as one frame. Based on the results from the atmospheric correction study over Hawaii and the findings over the deforming area, we applied the ECMWF correction or each of the frames we processed in the Yellowstone area. Figure 9 shows different atmospheric corrections applied for the 49-142 frame.
We used precise orbits since they are available to process past data. We show cumulative displacements between May 2017 and August 2018 for all frames, as some data were not available between January 2017 and April 2017 and between September 2018 and December 2018 ( Figure 10).   Table 1) obtained from four different acquisition frames of Sentinel-1 data. Red areas show uplift, which is apparently coinciding with steep valleys. Blue areas represent subsidence; black circles identify consistent subsidence on all four images.
These initial results show that some residual apparent uplift (red) is present in many cases, which is most likely due to a DEM or atmospheric correction error as they clearly coincide with steep topography gradients in the region. This apparent uplift marked with red is most prominent in the Path 49, Frame 142 image (top right panel on Figure 10). Subsidence in these images is marked with blue and we see a consistent signal in the southwest corner of the Yellowstone National Park on all four cumulative displacement images (black circles on Figure 10). This is the region where the Yellowstone volcano is most active, marked with geysers and hot springs and gas emission activities.
We collected and analysed the available continuous GNSS data and compared it to displacements derived from the DInSAR analysis. Figures 11-14 show comparison between (1) the displacement derived from GNSS time series in the LOS direction of the DInSAR acquisition and (2) the LOS displacement extracted from the DInSAR solution at the location corresponding to the GNSS station. In this comparison, we include all available DInSAR time series between January 2017 and December 2018. Some stations show better fit (MAWY, P720, or OFW2) than others (WLWY or P710). It is difficult to gain a general idea of the quality of the fit between GNSS and DInSAR displacement, because these permanent GNSS stations are present only at the area of the Yellowstone National Park, and most of the DInSAR frame do not have any ground truth data to compare with.

Discussion
We have compared the effects of applying two different atmospheric corrections, the ECMWF and GACOS corrections. The ECMWF correction is based on weather forecast and meteorological modelling only, while the GACOS also includes GNSS measurements. We found the two corrections similar in both cumulative displacement results and when we compared time series of individual pixels to GNSS data. An advantage of the ECMWF correction is that it is possible to extract the atmospheric information automatically through the processing, while the GACOS correction requires that we request the data and wait times can be up to 24 h. For the GeoSciFramework project purposes to provide early warning of changes in volcanic activities in the Yellowstone volcanic area, the ECMWF correction and its automated setup is better suited to the current computational structure. We will revisit this issue when GACOS becomes available on a similar automated basis.
Analysis at Hawaii for the time period January 2017 through June 2018 provides an assessment of the potential effect of using atmospheric corrections on hazard evaluation. Figure 4 shows the estimated change in height associated with atmospheric corrections in Hawaii, and the differences between the two methods. Note that the difference in estimated displacements between the ECMWF and no correction and the ECMWF and GACOS corrections range from 2.5 to 5 cm in those locations associated with either significant topography or magma motion. In the case of the former, this could result in overestimating potential hazard, while in the latter, it could mean inaccurate estimation of the size and location of the magma source, resulting in underestimation of the potential hazard. Future work will include assessment of these effects for a wider variety of hazard types and associated models.
We found that using real-time orbits instead of precise ones does not significantly affect the results of cumulative displacement over the Kilauea case study. There are pixel-to-pixel differences, but they do not change the overall results and interpretation of the DInSAR analysis. Therefore, it is possible to use these measurements in near-real time as they are acquired and processed, which is essential for early-warning purposes, as in the case of monitoring the Yellowstone volcanic area. As precise orbits are being continuously calculated and updated, the DInSAR displacement results can be refined for past acquisitions and real-time orbits employed for the most recent ones.
The important conclusion here is that the use of real-time orbits will not significantly impact the estimate of surface displacements (Figures 7 and 8) and, as such, not adversely affect the associated hazard estimates or evaluation, regardless of the hazard type.
Real-time orbits are available within two hours of the DInSAR image acquisition, and we expect to have updated ground motion and displacement results approximately eight hours after each new SAR image is available over the Yellowstone area. This estimation includes the time it takes to pull the new image data, wait the two hours for the real-time orbit to become available, complete the processing through GMT5SAR software, and apply the final GIAnT time series inversion. The new image would be interfered against only the most recent acquisitions contained within the 50-day threshold. This suggests that 6-12 new interferograms would be processed each time and then added to the continuously updating model of the Yellowstone volcanic area, together with other remote sensing observations from the GNSS network and ground-based borehole strain and seismic measurements. Through machine learning techniques, it is possible to ingest a large amount of data and set up a warning system that will provide alerts to any change in seismic activity or surface deformation.
The initial cumulative displacement calculations ( Figure 10) will contribute towards the baseline model that will be set up over the course of the GeoSciFramework project. Future work in the GeoSciFramework project includes implementing a faster processing chain described by Zebker [35], merging DInSAR observations from different acquisition frames [65] and optimizing integration of the different types of observations with their varying spatial and temporal acquisition rates.
It would be beneficial to install some GNSS stations across these frames outside of the Yellowstone National Park area, and acquire further ground truth data. This could contribute to the analysis of the displacements and could provide explanation to some of the apparent uplift that correlates with the topography.