Integration of DInSAR Time Series and GNSS Data for Continuous Volcanic Deformation Monitoring and Eruption Early Warning Applications

: With approximately 800 million people globally living within 100 km of a volcano, it is essential that we build a reliable observation system capable of delivering early warnings to potentially impacted nearby populations. Global Navigation Satellite System (GNSS) and satellite Synthetic Aperture Radar (SAR) document comprehensive ground motions or ruptures near, and at, the Earth’s surface and may be used to detect and analyze natural hazard phenomena. These datasets may also be combined to improve the accuracy of deformation results. Here, we prepare a differential interferometric SAR (DInSAR) time series and integrate it with GNSS data to create a fused dataset with enhanced accuracy of 3D ground motions over Hawaii island from November 2015 to April 2021. We present a comparison of the raw datasets against the fused time series and give a detailed account of observed ground deformation leading to the May 2018 and December 2020 volcanic eruptions. Our results provide important new estimates of the spatial and temporal dynamics of the 2018 Kilauea volcanic eruption. The methodology presented here can be easily repeated over any region of interest where an SAR scene overlaps with GNSS data. The results will contribute to diverse geophysical studies, including but not limited to the classiﬁcation of precursory movements leading to major eruptions and the advancement of early warning systems.


Introduction
Volcanic eruptions, earthquakes, and tsunamis occur over numerous spatial and temporal scales. Although these phenomena are often studied individually, there is frequently interconnectivity between disaster types. For example, concentrated swarms of earthquakes, elevated readings of gas emission, and increased ground motion over volcanic regions may indicate an impending eruption [1][2][3][4][5][6][7][8]. Most active volcanoes around the world are monitored using geodetic data sets such as Synthetic Aperture Radar (SAR) and Global Navigation Satellite System (GNSS) data, in conjunction with other ground-based instruments, with the goal of providing early warning for major eruptions and reducing risk to nearby populations or infrastructure [9]. While several studies have attempted to forecast or model potential volcano hazards using remote sensing techniques [2,[10][11][12][13][14][15][16], there is currently no single framework in place that simultaneously consolidates geodetic data from multiple sensors, freely provides scientists with near real-time continuous time series products and is capable of distinguishing and broadcasting geophysical events.
The GeoScience CyberInfrastructure Framework (GeoSCIFramework or GSF) project aims to improve intermediate-to-short term forecasts of catastrophic natural hazard events, allowing researchers to instantly detect phenomena and reveal more suppressed, longterm motions of Earth's surface at unprecedented spatial and temporal resolutions. These goals will be accomplished by applying big data analytics and training machine learning algorithms to recognize patterns across various data signals during noteworthy events. When complete, the system will be capable of processing and delivering large streams of near real-time data from a mix of Differential Interferometric SAR (DInSAR) imagery, GNSS, and other geodetic-related sensors, as well as seismic, gas emission, and thermal data.
DInSAR quantifies line-of-sight (LOS) ground deformation with mm-cm precision, and GNSS data delivers precise point positioning and timing data to determine exact location and deformation measurements, also with mm-cm scale precision. Furthermore, DInSAR processing can be combined with GNSS data to obtain 3D ground surface motions [17][18][19][20]. Together, these time series produce high resolution, sub-centimeter precision measurements of ground deformation over large swaths of Earth's surface with dense spatiotemporal coverage, which provides scientists with a greater understanding of crustal or shallow subsurface dynamics over volcanic regions.
We focus on generating an automated DInSAR time series processing routine that is integrated with GNSS data into a unified deformation field to provide more constrained deformation rates and vector measurements related to volcanic activity. We process Sentinel-1A/B SAR data into time series over Hawaii from November 2015 to April 2021 and integrate those results with GNSS data at various station positions ( Figure 1). The DInSAR + GNSS integrated time series can be used to describe the full extent of ground motions through time with decreased uncertainty in three directions of motion (east-west, north-south, and up-down). We present the unified DInSAR and GNSS time series and compare them to the original datasets.  [21]. The yellow box shows the cropped outline of each interferogram used when generating time series. Grey circles and colored diamonds indicate GNSS station locations, which were used to create kriging interpolated GNSS maps in Section 2.3. Twenty-four-hour final solution GNSS time series data, from stations listed in Table 1 and aligned to the local, fixed, Pacific Plate reference frame were obtained through the Nevada Geodetic Laboratory (NGL), University of Nevada Reno (http://geodesy.unr.edu/ (accessed on 1 November 2021)). Stations are maintained by the USGS HVO [22,23], and data are archived and distributed by the UNAVCO GAGE facility. We take a closer look at the time series over the diamond GNSS locations in Section 3.2, where the blue diamond is the NUPM GNSS station, the purple diamond corresponds to the CRIM station, orange represents the MKEA station, and green is the BLBP station. The blue triangle shows the location of the Mauna Loa volcano summit, and the yellow polygon indicates where the East Rift Zone is located. Background image taken from Google Earth/Data SIO, NOAA, U.S. Navy, NGA, GEBCO/Data LDEO-Columbia, NSF, NOAA/ Imagery Date: 13 December 2015.
Our region of interest is over the Big Island of Hawaii. Unlike most volcanic systems, Hawaii experiences frequent eruptive activity and provides an opportunity for scientists to record extensive observations over multiple events. Kilauea volcano has erupted 34 times since 1952 [24]. Two volcanic eruptions are captured within our time series between November 2015 and April 2021. The first eruption occurred in May 2018, when the lava lake within the Halema'uma'u crater (collocated with CRIM station, Figure 1) drained following an intrusion into, and subsequent eruption from, the Pu Remote Sens. 2022, 14, x FOR PEER REVIEW 4 of 33 following an intrusion into, and subsequent eruption from, the Puʻu ʻŌʻō crater (collocated with NUPM station, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow polygon, Figure 1). Sudden changes in lava lake levels, increased micro-seismicity around Kilauea's summit, and deformation along the NE strike of the ERZ were all precursory indications reported by the Hawaiian Volcano Observatory (HVO) that an impending eruption might occur [10,25]. According to [26], an obstruction in the magma plumbing system at Puʻu ʻŌʻō volcano caused widespread pressurization in the volcano, driving magma into the lower, eastern flank. Activity over the ERZ decelerated by September 2018 and remained quiet until December 2020, when a summit eruption that continued through May 2021 refilled the lava lake within the Halema'uma'u crater [24]. Although we observe this activity within the separate datasets, DInSAR and GNSS, as shown in Section 2, by integrating the two geodetic datasets together, we recover a 3-D spatial map, instead of 1-D LOS motion, with the resolution of the DInSAR images, that includes time series for motion in the east-west, north-south, and up direction, at each location. In addition, there is a significant improvement in the accuracy of each component of motion (east, north, up) relative to either the DInSAR or GNSS data alone, as presented in Sections 3 and 4.

Materials and Methods
The automated DInSAR processing routine was separated into three different components. The first section (Stage 1) was built based on GMTSAR source code [27,28] to process single look complex (SLC) DInSAR satellite imagery into interferograms using the small baseline subset (SBAS) method [29] performed in parallel. The second phase (Stage 2) applies the New SBAS (NSBAS) inversion [30] method to the GMTSAR interferograms of Stage 1 and generates the DInSAR time series and the cumulative deformation map. Finally, the third component (Stage 3) produces integrated 3D displacements using the LOS deformation from Stage 2, the geometry of the SAR acquisition, and precise, 3D vector positioning measurements [12,[18][19][20]. The routine can process the fused data at a single pixel, which is delivered as a plotted time series, or as an interpolated displacement map over a larger region, created from an array of available GNSS stations within the extent of the SAR scene.

Data
For this study, 250 descending Sentinel-1A/B SLC images were acquired between November 2015 to April 2021 along Path 87 Scene 526 ( Figure 1; Supplementary Material List S1) through the Alaskan Satellite Facility (ASF) Vertex portal [21]. Twenty-four-hour final solution GNSS data was managed by the USGS HVO and archived by the UNAVCO GAGE facility; processed time series were generated by and distributed through NGL [22,23]. Data from 48 GNSS stations over Hawaii (Table 1, Figure 1) were obtained over the same period of time, decimated to match the sampling rate of the InSAR time series, and used to create the interpolated map for this study. We used the ordinary kriging interpolation algorithm [31][32][33] supported by an exponential distribution model to construct standard variograms from the 48 GNSS station data.
While GNSS is known for its high precision in the horizontal directions (east and north), estimates of vertical motion have a larger uncertainty [34,35]. On the other hand, with an incidence angle range of 18.3° to 46.8°, DInSAR sensors are most sensitive to vertical displacements and can help to improve ground velocity estimates in the up direction [20,35]. Our integrated results provide a better representation, and therefore, a better understanding of the volcanic deformation pattern and subsurface-surface behavior through time. The workflow presented here can be easily repeated or applied to other locations where GNSS data falls within an SAR scene footprint. following an intrusion into, and subsequent eruption from, the Puʻu ʻŌʻō crater (collocated with NUPM station, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow polygon, Figure 1). Sudden changes in lava lake levels, increased micro-seismicity around Kilauea's summit, and deformation along the NE strike of the ERZ were all precursory indications reported by the Hawaiian Volcano Observatory (HVO) that an impending eruption might occur [10,25]. According to [26], an obstruction in the magma plumbing system at Puʻu ʻŌʻō volcano caused widespread pressurization in the volcano, driving magma into the lower, eastern flank. Activity over the ERZ decelerated by September 2018 and remained quiet until December 2020, when a summit eruption that continued through May 2021 refilled the lava lake within the Halema'uma'u crater [24]. Although we observe this activity within the separate datasets, DInSAR and GNSS, as shown in Section 2, by integrating the two geodetic datasets together, we recover a 3-D spatial map, instead of 1-D LOS motion, with the resolution of the DInSAR images, that includes time series for motion in the east-west, north-south, and up direction, at each location. In addition, there is a significant improvement in the accuracy of each component of motion (east, north, up) relative to either the DInSAR or GNSS data alone, as presented in Sections 3 and 4.

Materials and Methods
The automated DInSAR processing routine was separated into three different components. The first section (Stage 1) was built based on GMTSAR source code [27,28] to process single look complex (SLC) DInSAR satellite imagery into interferograms using the small baseline subset (SBAS) method [29] performed in parallel. The second phase (Stage 2) applies the New SBAS (NSBAS) inversion [30] method to the GMTSAR interferograms of Stage 1 and generates the DInSAR time series and the cumulative deformation map. Finally, the third component (Stage 3) produces integrated 3D displacements using the LOS deformation from Stage 2, the geometry of the SAR acquisition, and precise, 3D vector positioning measurements [12,[18][19][20]. The routine can process the fused data at a single pixel, which is delivered as a plotted time series, or as an interpolated displacement map over a larger region, created from an array of available GNSS stations within the extent of the SAR scene.

Data
For this study, 250 descending Sentinel-1A/B SLC images were acquired between November 2015 to April 2021 along Path 87 Scene 526 ( Figure 1; Supplementary Material List S1) through the Alaskan Satellite Facility (ASF) Vertex portal [21]. Twenty-four-hour final solution GNSS data was managed by the USGS HVO and archived by the UNAVCO GAGE facility; processed time series were generated by and distributed through NGL [22,23]. Data from 48 GNSS stations over Hawaii (Table 1, Figure 1) were obtained over the same period of time, decimated to match the sampling rate of the InSAR time series, and used to create the interpolated map for this study. We used the ordinary kriging interpolation algorithm [31][32][33] supported by an exponential distribution model to construct standard variograms from the 48 GNSS station data.
While GNSS is known for its high precision in the horizontal directions (east and north), estimates of vertical motion have a larger uncertainty [34,35]. On the other hand, with an incidence angle range of 18.3° to 46.8°, DInSAR sensors are most sensitive to vertical displacements and can help to improve ground velocity estimates in the up direction [20,35]. Our integrated results provide a better representation, and therefore, a better understanding of the volcanic deformation pattern and subsurface-surface behavior through time. The workflow presented here can be easily repeated or applied to other locations where GNSS data falls within an SAR scene footprint. following an intrusion into, and subsequent eruption from, the Puʻu ʻŌʻō crater (collocated with NUPM station, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow polygon, Figure 1). Sudden changes in lava lake levels, increased micro-seismicity around Kilauea's summit, and deformation along the NE strike of the ERZ were all precursory indications reported by the Hawaiian Volcano Observatory (HVO) that an impending eruption might occur [10,25]. According to [26], an obstruction in the magma plumbing system at Puʻu ʻŌʻō volcano caused widespread pressurization in the volcano, driving magma into the lower, eastern flank. Activity over the ERZ decelerated by September 2018 and remained quiet until December 2020, when a summit eruption that continued through May 2021 refilled the lava lake within the Halema'uma'u crater [24]. Although we observe this activity within the separate datasets, DInSAR and GNSS, as shown in Section 2, by integrating the two geodetic datasets together, we recover a 3-D spatial map, instead of 1-D LOS motion, with the resolution of the DInSAR images, that includes time series for motion in the east-west, north-south, and up direction, at each location. In addition, there is a significant improvement in the accuracy of each component of motion (east, north, up) relative to either the DInSAR or GNSS data alone, as presented in Sections 3 and 4.

Materials and Methods
The automated DInSAR processing routine was separated into three different components. The first section (Stage 1) was built based on GMTSAR source code [27,28] to process single look complex (SLC) DInSAR satellite imagery into interferograms using the small baseline subset (SBAS) method [29] performed in parallel. The second phase (Stage 2) applies the New SBAS (NSBAS) inversion [30] method to the GMTSAR interferograms of Stage 1 and generates the DInSAR time series and the cumulative deformation map. Finally, the third component (Stage 3) produces integrated 3D displacements using the LOS deformation from Stage 2, the geometry of the SAR acquisition, and precise, 3D vector positioning measurements [12,[18][19][20]. The routine can process the fused data at a single pixel, which is delivered as a plotted time series, or as an interpolated displacement map over a larger region, created from an array of available GNSS stations within the extent of the SAR scene.

Data
For this study, 250 descending Sentinel-1A/B SLC images were acquired between November 2015 to April 2021 along Path 87 Scene 526 ( Figure 1; Supplementary Material List S1) through the Alaskan Satellite Facility (ASF) Vertex portal [21]. Twenty-four-hour final solution GNSS data was managed by the USGS HVO and archived by the UNAVCO GAGE facility; processed time series were generated by and distributed through NGL [22,23]. Data from 48 GNSS stations over Hawaii (Table 1, Figure 1) were obtained over the same period of time, decimated to match the sampling rate of the InSAR time series, and used to create the interpolated map for this study. We used the ordinary kriging interpolation algorithm [31][32][33] supported by an exponential distribution model to construct standard variograms from the 48 GNSS station data.
While GNSS is known for its high precision in the horizontal directions (east and north), estimates of vertical motion have a larger uncertainty [34,35]. On the other hand, with an incidence angle range of 18.3° to 46.8°, DInSAR sensors are most sensitive to vertical displacements and can help to improve ground velocity estimates in the up direction [20,35]. Our integrated results provide a better representation, and therefore, a better understanding of the volcanic deformation pattern and subsurface-surface behavior through time. The workflow presented here can be easily repeated or applied to other locations where GNSS data falls within an SAR scene footprint. ō crater (collocated with NUPM station, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow polygon, Figure 1). Sudden changes in lava lake levels, increased micro-seismicity around Kilauea's summit, and deformation along the NE strike of the ERZ were all precursory indications reported by the Hawaiian Volcano Observatory (HVO) that an impending eruption might occur [10,25]. According to [26], an obstruction in the magma plumbing system at Pu 4 of 33 , and subsequent eruption from, the Puʻu ʻŌʻō crater (collo- Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow changes in lava lake levels, increased micro-seismicity around rmation along the NE strike of the ERZ were all precursory Hawaiian Volcano Observatory (HVO) that an impending ]. According to [26], an obstruction in the magma plumbing o caused widespread pressurization in the volcano, driving rn flank. Activity over the ERZ decelerated by September 2018 cember 2020, when a summit eruption that continued through ke within the Halema'uma'u crater [24]. Although we observe rate datasets, DInSAR and GNSS, as shown in Section 2, by datasets together, we recover a 3-D spatial map, instead of 1olution of the DInSAR images, that includes time series for th-south, and up direction, at each location. In addition, there t in the accuracy of each component of motion (east, north, SAR or GNSS data alone, as presented in Sections 3 and 4. R processing routine was separated into three different comtage 1) was built based on GMTSAR source code [27,28] to (SLC) DInSAR satellite imagery into interferograms using the ) method [29] performed in parallel. The second phase (Stage SBAS) inversion [30] method to the GMTSAR interferograms e DInSAR time series and the cumulative deformation map. t (Stage 3) produces integrated 3D displacements using the e 2, the geometry of the SAR acquisition, and precise, 3D vects [12,[18][19][20]. The routine can process the fused data at a sinas a plotted time series, or as an interpolated displacement reated from an array of available GNSS stations within the ending Sentinel-1A/B SLC images were acquired between Nolong Path 87 Scene 526 ( Figure 1; Supplementary Material List ellite Facility (ASF) Vertex portal [21]. Twenty-four-hour final anaged by the USGS HVO and archived by the UNAVCO me series were generated by and distributed through NGL stations over Hawaii (Table 1, Figure 1) were obtained over cimated to match the sampling rate of the InSAR time series, polated map for this study. We used the ordinary kriging in-] supported by an exponential distribution model to construct he 48 GNSS station data.
for its high precision in the horizontal directions (east and motion have a larger uncertainty [34,35]. On the other hand, ge of 18.3° to 46.8°, DInSAR sensors are most sensitive to verhelp to improve ground velocity estimates in the up direction to, and subsequent eruption from, the Puʻu ʻŌʻō crater (collon, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow n changes in lava lake levels, increased micro-seismicity around eformation along the NE strike of the ERZ were all precursory the Hawaiian Volcano Observatory (HVO) that an impending ,25]. According to [26], an obstruction in the magma plumbing ano caused widespread pressurization in the volcano, driving stern flank. Activity over the ERZ decelerated by September 2018 December 2020, when a summit eruption that continued through lake within the Halema'uma'u crater [24]. Although we observe eparate datasets, DInSAR and GNSS, as shown in Section 2, by tic datasets together, we recover a 3-D spatial map, instead of 1resolution of the DInSAR images, that includes time series for orth-south, and up direction, at each location. In addition, there ent in the accuracy of each component of motion (east, north, DInSAR or GNSS data alone, as presented in Sections 3 and 4.
s AR processing routine was separated into three different com-(Stage 1) was built based on GMTSAR source code [27,28] to lex (SLC) DInSAR satellite imagery into interferograms using the AS) method [29] performed in parallel. The second phase (Stage (NSBAS) inversion [30] method to the GMTSAR interferograms the DInSAR time series and the cumulative deformation map. nent (Stage 3) produces integrated 3D displacements using the age 2, the geometry of the SAR acquisition, and precise, 3D vecents [12,[18][19][20]. The routine can process the fused data at a sinred as a plotted time series, or as an interpolated displacement , created from an array of available GNSS stations within the escending Sentinel-1A/B SLC images were acquired between No-1 along Path 87 Scene 526 ( Figure 1; Supplementary Material List atellite Facility (ASF) Vertex portal [21]. Twenty-four-hour final managed by the USGS HVO and archived by the UNAVCO time series were generated by and distributed through NGL SS stations over Hawaii (Table 1, Figure 1) were obtained over decimated to match the sampling rate of the InSAR time series, terpolated map for this study. We used the ordinary kriging in--33] supported by an exponential distribution model to construct the 48 GNSS station data. wn for its high precision in the horizontal directions (east and cal motion have a larger uncertainty [34,35]. On the other hand, ange of 18.3° to 46.8°, DInSAR sensors are most sensitive to veran help to improve ground velocity estimates in the up direction into, and subsequent eruption from, the Puʻu ʻŌʻō crater (collotion, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow den changes in lava lake levels, increased micro-seismicity around deformation along the NE strike of the ERZ were all precursory y the Hawaiian Volcano Observatory (HVO) that an impending [10,25]. According to [26], an obstruction in the magma plumbing olcano caused widespread pressurization in the volcano, driving eastern flank. Activity over the ERZ decelerated by September 2018 til December 2020, when a summit eruption that continued through va lake within the Halema'uma'u crater [24]. Although we observe separate datasets, DInSAR and GNSS, as shown in Section 2, by detic datasets together, we recover a 3-D spatial map, instead of 1he resolution of the DInSAR images, that includes time series for t, north-south, and up direction, at each location. In addition, there ement in the accuracy of each component of motion (east, north, e DInSAR or GNSS data alone, as presented in Sections 3 and 4.

ods
InSAR processing routine was separated into three different comion (Stage 1) was built based on GMTSAR source code [27,28] to plex (SLC) DInSAR satellite imagery into interferograms using the SBAS) method [29] performed in parallel. The second phase (Stage S (NSBAS) inversion [30] method to the GMTSAR interferograms tes the DInSAR time series and the cumulative deformation map. ponent (Stage 3) produces integrated 3D displacements using the Stage 2, the geometry of the SAR acquisition, and precise, 3D vecements [12,[18][19][20]. The routine can process the fused data at a sinivered as a plotted time series, or as an interpolated displacement ion, created from an array of available GNSS stations within the e.
descending Sentinel-1A/B SLC images were acquired between No-021 along Path 87 Scene 526 ( Figure 1; Supplementary Material List n Satellite Facility (ASF) Vertex portal [21]. Twenty-four-hour final as managed by the USGS HVO and archived by the UNAVCO sed time series were generated by and distributed through NGL NSS stations over Hawaii (Table 1, Figure 1) were obtained over e, decimated to match the sampling rate of the InSAR time series, interpolated map for this study. We used the ordinary kriging in-31-33] supported by an exponential distribution model to construct rom the 48 GNSS station data. nown for its high precision in the horizontal directions (east and rtical motion have a larger uncertainty [34,35]. On the other hand, e range of 18.3° to 46.8°, DInSAR sensors are most sensitive to verd can help to improve ground velocity estimates in the up direction ō volcano caused widespread pressurization in the volcano, driving magma into the lower, eastern flank. Activity over the ERZ decelerated by September 2018 and remained quiet until December 2020, when a summit eruption that continued through May 2021 refilled the lava lake within the Halema'uma'u crater [24]. Although we observe this activity within the separate datasets, DInSAR and GNSS, as shown in Section 2, by integrating the two geodetic datasets together, we recover a 3-D spatial map, instead of 1-D LOS motion, with the resolution of the DInSAR images, that includes time series for motion in the east-west, north-south, and up direction, at each location. In addition, there is a significant improvement in the accuracy of each component of motion (east, north, up) relative to either the DInSAR or GNSS data alone, as presented in Sections 3 and 4.

Materials and Methods
The automated DInSAR processing routine was separated into three different components. The first section (Stage 1) was built based on GMTSAR source code [27,28] to process single look complex (SLC) DInSAR satellite imagery into interferograms using the small baseline subset (SBAS) method [29] performed in parallel. The second phase (Stage 2) applies the New SBAS (NSBAS) inversion [30] method to the GMTSAR interferograms of Stage 1 and generates the DInSAR time series and the cumulative deformation map. Finally, the third component (Stage 3) produces integrated 3D displacements using the LOS deformation from Stage 2, the geometry of the SAR acquisition, and precise, 3D vector positioning measurements [12,[18][19][20]. The routine can process the fused data at a single pixel, which is delivered as a plotted time series, or as an interpolated displacement map over a larger region, created from an array of available GNSS stations within the extent of the SAR scene.

Data
For this study, 250 descending Sentinel-1A/B SLC images were acquired between November 2015 to April 2021 along Path 87 Scene 526 (Figure 1; Supplementary Material List S1) through the Alaskan Satellite Facility (ASF) Vertex portal [21]. Twenty-four-hour final solution GNSS data was managed by the USGS HVO and archived by the UN-AVCO GAGE facility; processed time series were generated by and distributed through NGL [22,23]. Data from 48 GNSS stations over Hawaii (Table 1, Figure 1) were obtained over the same period of time, decimated to match the sampling rate of the InSAR time series, and used to create the interpolated map for this study. We used the ordinary kriging interpolation algorithm [31][32][33] supported by an exponential distribution model to construct standard variograms from the 48 GNSS station data.
While GNSS is known for its high precision in the horizontal directions (east and north), estimates of vertical motion have a larger uncertainty [34,35]. On the other hand, with an incidence angle range of 18.3 • to 46.8 • , DInSAR sensors are most sensitive to vertical displacements and can help to improve ground velocity estimates in the up direction [20,35]. Our integrated results provide a better representation, and therefore, a better understanding of the volcanic deformation pattern and subsurface-surface behavior through time. The workflow presented here can be easily repeated or applied to other locations where GNSS data falls within an SAR scene footprint. Table 1. GNSSS stations over Hawaii Island used in this study for time series integration with DInSAR data. Bolded station corresponds to colored diamonds in Figure 1. GNSS data were obtained through the NGL, and stations were maintained by the USGS HVO [22,23].

Building Interferograms
DInSAR processing is challenging and complex in that several improvements, including orbital, topographic, and atmospheric corrections, must be applied to isolate the deformation signal. As we consider an automated system that will be streaming continuous Sentinel-1A/B data over Hawaii Island, we must decide how these corrections are applied and systematized within the processing routine. This first component requires several input parameters and files to run, including the SAR images, the corresponding orbital correction files, a Digital Elevation Model (DEM), and the number of allowed days between potential interferometric image pairs.
For the time series presented in this paper, we used only precise orbital corrections on each Sentinel-1 A/B SLC image. These precise files, which were also obtained through the ASF Vertex portal [21], are not available to users until two weeks after the SAR image is acquired by the Sentinel satellite. Thus, for the real-time automated system, and due to the Sentinel-1A/B satellites' repeat period of 6-12 days, the most recent SAR image acquisitions over the region of interest were first processed using real-time orbital corrections and then updated with the precise corrections two weeks later. We thoroughly review the effects of real-time vs. precise orbital corrections in [9] and determine that the use of real-time orbits is sufficient for early-warning applications. For the topographic correction, we used a 30 m resolution SRTM DEM that completely covers the SAR image footprint.
After the orbital and topographic corrections are applied, the complex interferogram is formed. We allowed the algorithm to pair any two images that were obtained within 35 days of each other. From the 250 descending SLC images acquired between November 2015 and April 2021, 671 interferometric image pairs were successfully generated (Supplementary Material List S2). For each pair, interferometric products of phase, coherence, phase gradient, and LOS displacement were constructed in both radar and geographical coordinates ( Figure 2). The code utilizes the Snaphu phase unwrapping method [36] and geocodes all interferometric output products. These outputs are formatted as GMT-compatible grids, which easily convert to GeoTIFF or other GDAL-compatible raster drivers.

DInSAR Time Series Generation
The second component of the automated system uses the GIAnT software [37] to construct LOS displacement maps and time series. In GIAnT, unwrapped interferograms are converted to units of millimeters before any processing takes place. The software assembles the input data, the unwrapped interferograms, coherence maps, masks, and metadata, into binary files and a composite HDF5 file that is compatible with GIAnT processing.
The user can choose between several atmospheric model corrections, including the European Centre for Medium-Range Weather Forecasts (ECMWF) model, which is built into the software so that the atmospheric files are pulled and applied automatedly. Currently, ECMWF is in the midst of migrating its data to an updated server, rendering the files temporarily unavailable. Therefore, for this study, the Generic Atmospheric Correction Online Service (GACOS) atmospheric model correction was applied. The GACOS model integrates the ECMWF atmospheric files with continuous GNSS tropospheric delay estimates, which currently must be manually ordered through the online portal [38][39][40][41]. The goal is to incorporate it into automated processing in the future.
The user can also choose between different time series inversion methods, available from GIAnT, to estimate filtered time series. The SBAS and NSBAS techniques, which are implemented here, are applicable when the differential interferograms have a small spatial baseline, a characteristic of Sentinel-1A/B data [29,30]. In particular, the NSBAS technique also estimates DEM errors and compensates for pixels that have missing observations. This inversion method estimates the LOS phase change of each pixel independently using a linear system.
For Stage 1, the automated scripts ran using high-performance computing nodes on the SUMMIT supercomputer located at the University of Colorado, Boulder, for approximately eight hours to produce 694 interferograms between November 2015 and April 2021. The NSBAS technique applies inversion on the 694 input interferograms to generate a time series with a cumulative displacement measurement for every date, or time-slice, included within the data. In this case, the time series produced 275 time-slices (See Figure 3). The atmospheric correction and inversion process in Stage 2 took approximately two hours to complete. Once the historical data has been processed, however, adding a new acquisition to the time series takes on the order of a couple of hours to fully acquire and process the necessary data, then output all the updated DInSAR-related products and time series.

Integration of Geodetic Datasets
We integrate our DInSAR time series with GNSS data from 48 stations over Hawaii to produce 3D high-resolution cumulative displacement maps with corresponding errors. The 24 h final GNSS solutions are provided in three components (east, north, up) and aligned to the local, fixed Pacific plate reference frame to minimize linear trends of the tectonic plate motion. It is worth noting that NGL also provides solutions aligned to the International GNSS Service-14 (IGS14) reference frame, which is based on the International Terrestrial Reference Frame (ITRF) and holds a no-net-rotation (NNR) [42]. We present a comparison of the results from both reference frames in the Supplementary Material (Figures S1-S3). First, we temporally subset the GNSS data from 48 stations to extract the data corresponding to the interval of the DInSAR time series. For each time step, we used the ordinary kriging algorithm and the 3D displacements from the 48 stations to compute a variogram. The spatial covariance of the GNSS data determines the structure of the variogram, and the weights were calculated from the data using an exponential distribution model to interpolate undefined points across the spatial field [31,32]. The kriging algorithm outputs three cumulative displacement maps (east, north, and up) with the same discretization and geocoding as the DInSAR output from Section 2.2 (Figure 4), [12,18,19,[43][44][45]. To combine all the geodetic displacement data into a single time series product, we expand on the DInSAR+GNSS velocity integration method developed by [12]. The method uses a Bayesian statistical modeling approach in conjunction with Markov random field (MRF) theory to combine datasets with varying spatio-temporal extents. Integration of DIn-SAR and GNSS high-resolution, 3D surface displacement measurements was achieved by minimization of the energy function related to the corresponding Gibbs random field (GRF) distribution, in which the joint probability density of the variables are positive [12,18,19,33]. The energy function is as follows: where U b a is the likelihood energy, b is the observation with uncertainty σ, a is the unknown parameter, and N is the number of observations or pixels within the acquisition. The resulting adaptation of Equation (1) to our geodetic displacement datasets is: with coefficients: where σ is the standard deviation for the measurements, d LOS is the cumulative LOS displacement, S LOS Uncertainties associated with the 3D displacement maps and the plotted time series are automatically generated with the corresponding products. These are integrated using the same methodology as the observations and are a combination of the raw GNSS data uncertainties and the LOS DInSAR uncertainties (Equations S6-S8 in Supplementary Material). For the 3D displacement maps, an additional error is introduced when we interpolate data from the 48 GNSS stations into a variogram using the ordinary kriging algorithm. We present the error analysis for that integration in the Supplementary Material ( Figures S4 and S5).

Cumulative Deformation Maps
The automated processing system can provide the assimilated DInSAR+GNSS time series to users in multiple ways. Over large areas, the code produced 3D high-resolution cumulative displacement maps and the corresponding uncertainties for each individual date or time-step of the series. Figure 5 shows the evolution of 3D surface motions captured in our time series, revealing the long-term deformation response to the May 2018 and December 2020 eruptions. The pixel resolution of the fused products was 100 m, the same as the DInSAR results.
Slight inflation (~10 cm) over the Kilauea volcano (purple diamond, Figures 1 and 5) was observed in the fused DInSAR+GNSS results as early as November 2016, nearly 1.5 years before the eruption took place (Figure 5a). By August 2017 (Figure 5b), cumulative inflation over Kilauea reached +26 cm. This magma intrusion is further supported by the pattern observed over the same region in the north-component subplot, in which the southern flank of Kilauea moved approximately 15 cm further south, and the northern flank moved 8 cm further north.
Upward motion over Mauna Loa volcano (blue triangle, Figures 1 and 5) also began to increase in August 2017. By February 2018 (Figure 5c), this inflation pattern was recognizable over Mauna Loa, where the ground reached +16 cm of uplift and was obvious over the Kilauea volcano after reaching +30 cm of total uplift. Over the following months, inflation at Mauna Loa began to slow, while it continued to increase at a similar/opposite rate over Kilauea. This interaction between the two volcanic systems has been studied before [46][47][48][49] and provides further evidence of upper mantle links and magma transportation. On 30 April 2018, the Pu'u 'Ō'ō eruptive vent collapsed [48][49][50], and on May 4, a Mw 6.9 earthquake occurred on the south flank of Kilauea, initiating the 2018 volcanic eruption. Our integrated results from 5 May 2018 (Figure 5d) recorded the 43 cm southeast ground rupture of the earthquake. Next, the sudden evacuation of subsurface material resulted in a rapid −10 cm ground deflation at Kilauea's crater and −80 cm over the ERZ near Pu Remote Sens. 2022, 14, x FOR PEER REVIEW following an intrusion into, and subsequent eruption from, the Puʻu ʻŌʻō cra cated with NUPM station, Figure 1) and Kilauea's lower East Rift Zone (ERZ polygon, Figure 1). Sudden changes in lava lake levels, increased micro-seismic Kilauea's summit, and deformation along the NE strike of the ERZ were all p indications reported by the Hawaiian Volcano Observatory (HVO) that an i eruption might occur [10,25]. According to [26], an obstruction in the magma system at Puʻu ʻŌʻō volcano caused widespread pressurization in the volcan magma into the lower, eastern flank. Activity over the ERZ decelerated by Septe and remained quiet until December 2020, when a summit eruption that continue May 2021 refilled the lava lake within the Halema'uma'u crater [24]. Although w this activity within the separate datasets, DInSAR and GNSS, as shown in Se integrating the two geodetic datasets together, we recover a 3-D spatial map, in u 4 of 33 d subsequent eruption from, the Puʻu ʻŌʻō crater (colloure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow nges in lava lake levels, increased micro-seismicity around ation along the NE strike of the ERZ were all precursory awaiian Volcano Observatory (HVO) that an impending ccording to [26], an obstruction in the magma plumbing aused widespread pressurization in the volcano, driving flank. Activity over the ERZ decelerated by September 2018 ber 2020, when a summit eruption that continued through within the Halema'uma'u crater [24]. Although we observe te datasets, DInSAR and GNSS, as shown in Section 2, by Ō 4 of 33 and subsequent eruption from, the Puʻu ʻŌʻō crater (collo- Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow hanges in lava lake levels, increased micro-seismicity around rmation along the NE strike of the ERZ were all precursory Hawaiian Volcano Observatory (HVO) that an impending ]. According to [26], an obstruction in the magma plumbing o caused widespread pressurization in the volcano, driving rn flank. Activity over the ERZ decelerated by September 2018 cember 2020, when a summit eruption that continued through ke within the Halema'uma'u crater [24]. Although we observe rate datasets, DInSAR and GNSS, as shown in Section 2, by ō volcano (yellow polygon, Figures 1 and 5a). As this eruptive activity began over the ERZ, the ground at Mauna Loa inflated in the east-west directions from a continued subsurface magma injection. This divergent, horizontal pattern became more obvious by June 2018 (Figure 5e), after Mauna Loa experienced an additional +3 cm of uplift (a net total of 13 cm). A closer look at deformation over the Kilauea crater in June 2018 indicated additional, sudden ruptures that moved the localized ground 25 cm to the west and in the vertical direction, Kilauea crater subsided another −9 cm. This deformation was likely due to the accelerated down-drop of Kilauea's caldera combined with 62 summit collapse events between May 16 and August 2 [48,49]. HVO reported that some collapse events released energy equivalent to a Mw 4.7 to 5.4 earthquake; our results captured these details.
Another interesting feature includes the three zones of subsidence radially surrounding Kilauea's summit and the Halema'uma'u crater. These deflated regions are spatially separated by~120 degrees. The two zones to the north and the southwest both experienced between −35 and −45 cm subsidence. The southeast leg of the three-pronged sinking feature fed into the ERZ, which has subsided more than −1.17 m to date. No active surface lava was observed after 4 September 2018 [48,49], but the ground surface over the ERZ continued to deflate through December of that year (Figure 5f), reaching a cumulative minimum of −1.32 m.
Soon after, the cumulative subsidence from the ERZ began to shrink, indicating that during the recovery from the 2018 eruption, the subsurface chambers may have begun refilling with magma. In November 2019 (Figure 5g), the east-west inflation pattern over Mauna Loa became more prominent, and uplift reached a total of +21 cm. The north-south component revealed that the Kilauea crater also recovered 17 cm of motion to the north (against the dominant southeastern motion of the ERZ and the coastline just south of the region) and regained +26 cm in vertical surface height.
Over the next year leading to the December 2020 summit eruption, another +11 cm of upward motion occurred at Mauna Loa (net total of +32 cm), and +14 cm occurred along the ERZ (net total of −83 cm) (Figure 5h). Immediately following the summit eruption, the ground height at Kilauea summit measured a total of +20 cm, while the region surrounding the crater subsided once more, decreasing the total ground height by −30 cm. In mid-April 2021, when our time series ended ( Figure 5I), a cumulative total of +21 cm of uplift occurred at the Kilauea crater, and +46 cm of uplift had occurred at Mauna Loa in response to the 2018 eruption and 2020 summit eruptions. Kilauea's surrounding volcanic flank and the ERZ experienced −50 cm and −1.0 m of subsidence, respectively.
These products allowed us to observe the total cumulative deformation pattern in the east, north, and up directions over the entire island from November 2015 to April 2021 with improved accuracy and detail. The 3D cumulative displacement maps provided new information regarding the pre-, during-, and post-eruption phases of the Hawaiian volcanic system at unprecedented spatial scales and revealed surface effects from magma movement and seismic activity leading to two different types of eruptions. For example, Figure 6 converts the integrated results and raw GNSS data from 3D to 1D, LOS deformation to match the raw DInSAR results, for easy comparison of all final results across the three datasets. While all three datasets visually complement one another, the results from integrating the GNSS data with the DInSAR data, Figure 6A, showed a more constrained uplift pattern than the DInSAR or GNSS results alone. Figure 7 compares the integrated results with the kriged GNSS data, illustrating the additional information provided by incorporating the DInSAR data. Uncertainty estimates for each component of motion at each time step are provided and analyzed in the Supplementary Material (Figures S6-S10).

Plotted Time Series
In addition to generating the assimilated displacement maps, we present time series at four locations within the image where the GNSS stations overlap with a corresponding DInSAR pixel coordinate. These stations are strategically positioned across the island and include both active and relatively stable regions. Figures 8-11 present the DInSAR+GNSS integrated time series at CRIM, NUPM, MKEA, and BLBP GNSS stations, respectively. We compare each component of motion from the integrated dataset with the raw GNSS time series and LOS DInSAR time series. These graphs depict the ground movement of a single point at every time-step of the series and record the surface response of that specific location to the two volcanic eruptions. The selected GNSS stations monitor a mixture of the most and least active regions of Hawaii and span a large spatial extent of the island. Similar to the 3D displacement maps, these individual time series provide researchers with a comprehensive study of surface deformation through time. By comparing the integrated product against the original GNSS time series, we find that information in the horizontal direction is slightly improved and information in the vertical direction is significantly enhanced once fused with LOS DInSAR measurements. The plotted time series provide additional insight into the temporal behavior at any given location.
The CRIM station (Figure 8) is located at the summit of Hawaii's most active volcano, Kilauea, along the southern rim of the Halema'uma'u crater (purple diamond, Figure 1 following an intrusion into, and subsequent eruption from, the Puʻu ʻŌʻō crater (collocated with NUPM station, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow polygon, Figure 1). Sudden changes in lava lake levels, increased micro-seismicity around Kilauea's summit, and deformation along the NE strike of the ERZ were all precursory indications reported by the Hawaiian Volcano Observatory (HVO) that an impending eruption might occur [10,25]. According to [26], an obstruction in the magma plumbing system at Puʻu ʻŌʻō volcano caused widespread pressurization in the volcano, driving magma into the lower, eastern flank. Activity over the ERZ decelerated by September 2018 and remained quiet until December 2020, when a summit eruption that continued through May 2021 refilled the lava lake within the Halema'uma'u crater [24]. Although we observe this activity within the separate datasets, DInSAR and GNSS, as shown in Section 2, by integrating the two geodetic datasets together, we recover a 3-D spatial map, instead of 1-D LOS motion, with the resolution of the DInSAR images, that includes time series for motion in the east-west, north-south, and up direction, at each location. In addition, there is a significant improvement in the accuracy of each component of motion (east, north, up) relative to either the DInSAR or GNSS data alone, as presented in Sections 3 and 4.

Materials and Methods
The automated DInSAR processing routine was separated into three different components. The first section (Stage 1) was built based on GMTSAR source code [27,28] to process single look complex (SLC) DInSAR satellite imagery into interferograms using the small baseline subset (SBAS) method [29] performed in parallel. The second phase (Stage 2) applies the New SBAS (NSBAS) inversion [30] method to the GMTSAR interferograms of Stage 1 and generates the DInSAR time series and the cumulative deformation map. Finally, the third component (Stage 3) produces integrated 3D displacements using the LOS deformation from Stage 2, the geometry of the SAR acquisition, and precise, 3D vector positioning measurements [12,[18][19][20]. The routine can process the fused data at a single pixel, which is delivered as a plotted time series, or as an interpolated displacement map over a larger region, created from an array of available GNSS stations within the extent of the SAR scene. following an intrusion into, and subsequent eruption from, the Puʻu ʻŌʻō crater (collocated with NUPM station, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow polygon, Figure 1). Sudden changes in lava lake levels, increased micro-seismicity around Kilauea's summit, and deformation along the NE strike of the ERZ were all precursory indications reported by the Hawaiian Volcano Observatory (HVO) that an impending eruption might occur [10,25]. According to [26], an obstruction in the magma plumbing system at Puʻu ʻŌʻō volcano caused widespread pressurization in the volcano, driving magma into the lower, eastern flank. Activity over the ERZ decelerated by September 2018 and remained quiet until December 2020, when a summit eruption that continued through May 2021 refilled the lava lake within the Halema'uma'u crater [24]. Although we observe this activity within the separate datasets, DInSAR and GNSS, as shown in Section 2, by integrating the two geodetic datasets together, we recover a 3-D spatial map, instead of 1-D LOS motion, with the resolution of the DInSAR images, that includes time series for motion in the east-west, north-south, and up direction, at each location. In addition, there is a significant improvement in the accuracy of each component of motion (east, north, up) relative to either the DInSAR or GNSS data alone, as presented in Sections 3 and 4.

Materials and Methods
The automated DInSAR processing routine was separated into three different components. The first section (Stage 1) was built based on GMTSAR source code [27,28] to process single look complex (SLC) DInSAR satellite imagery into interferograms using the small baseline subset (SBAS) method [29] performed in parallel. The second phase (Stage 2) applies the New SBAS (NSBAS) inversion [30] method to the GMTSAR interferograms of Stage 1 and generates the DInSAR time series and the cumulative deformation map. Finally, the third component (Stage 3) produces integrated 3D displacements using the LOS deformation from Stage 2, the geometry of the SAR acquisition, and precise, 3D vector positioning measurements [12,[18][19][20]. The routine can process the fused data at a single pixel, which is delivered as a plotted time series, or as an interpolated displacement map over a larger region, created from an array of available GNSS stations within the extent of the SAR scene. following an intrusion into, and subsequent eruption from, the Puʻu ʻŌʻō crater (collocated with NUPM station, Figure 1) and Kilauea's lower East Rift Zone (ERZ) (Yellow polygon, Figure 1). Sudden changes in lava lake levels, increased micro-seismicity around Kilauea's summit, and deformation along the NE strike of the ERZ were all precursory indications reported by the Hawaiian Volcano Observatory (HVO) that an impending eruption might occur [10,25]. According to [26], an obstruction in the magma plumbing system at Puʻu ʻŌʻō volcano caused widespread pressurization in the volcano, driving magma into the lower, eastern flank. Activity over the ERZ decelerated by September 2018 and remained quiet until December 2020, when a summit eruption that continued through May 2021 refilled the lava lake within the Halema'uma'u crater [24]. Although we observe this activity within the separate datasets, DInSAR and GNSS, as shown in Section 2, by integrating the two geodetic datasets together, we recover a 3-D spatial map, instead of 1-D LOS motion, with the resolution of the DInSAR images, that includes time series for motion in the east-west, north-south, and up direction, at each location. In addition, there is a significant improvement in the accuracy of each component of motion (east, north, up) relative to either the DInSAR or GNSS data alone, as presented in Sections 3 and 4.

Materials and Methods
The automated DInSAR processing routine was separated into three different components. The first section (Stage 1) was built based on GMTSAR source code [27,28] to process single look complex (SLC) DInSAR satellite imagery into interferograms using the small baseline subset (SBAS) method [29] performed in parallel. The second phase (Stage 2) applies the New SBAS (NSBAS) inversion [30] method to the GMTSAR interferograms of Stage 1 and generates the DInSAR time series and the cumulative deformation map. Finally, the third component (Stage 3) produces integrated 3D displacements using the LOS deformation from Stage 2, the geometry of the SAR acquisition, and precise, 3D vector positioning measurements [12,[18][19][20]. The routine can process the fused data at a single pixel, which is delivered as a plotted time series, or as an interpolated displacement map over a larger region, created from an array of available GNSS stations within the extent of the SAR scene. ō volcano (blue diamond, Figure 1), which connects with the ERZ. We present plotted time series from this location to show the sudden displacement from the May 4 Mw 6.9 earthquake. The ground at NUPM stayed relatively stable between November 2015 and May 2018, leading up to the volcanic eruption, with~4.5 cm southeast motion and stable (+/−1 cm) fluctuations in the vertical direction. Suddenly, the signal broke 58.83 cm southeast and 28 cm downward. Volcanic activity took over the signal and continued over the next couple of months, adding 20.88 cm of horizontal motion to the southeast (net total of 79.71 cm) and another 38 cm of subsidence (net total of 66 cm) before entering a state of ground recovery. Leading to the December 2020 summit eruption, the ground at the NUPM station recovered~8.95 cm to the northwest and regained 17 cm of uplift. The summit eruption was slightly detectable at this location, having recorded 2 cm of motion to the northeast.
The MKEA station ( Figure 10) lies along the southeastern flank of the dormant volcano, Mauna Kea (orange diamond, Figure 1). Data from November 2015 to January 2020 at the MKEA station was associated with the highest uncertainty values throughout this study. In the horizontal directions, the integrated dataset agreed nicely with the GNSS data; however, in the vertical orientation, the data trended more with the DInSAR time series results between 2015 and 2019. In January 2020, HVO reported an instrumental dome replacement, after which the variability in the integrated signal tightened, yet still deviated slightly (~1 cm) from the raw GNSS data. This time series provides an important example of the GNSS dataset with the highest error used in this study and demonstrates how combing the GNSS data with DInSAR also keeps a system of checks and balances of the integrated system. Even with high uncertainty, we could distinguish a 2.5 cm trend southeastward during the 2018 eruption.    The scatter observed in the DInSAR data, and the integrated results in the up direction in Figure 10 are partially due to the fact that the GPS data is at a single, isolated point, while the DInSAR results average the value of deformation over a 100 × 100 m 2 pixel. The DInSAR values may also contain small amounts of residual tropospheric noise. Examination of the MKEA station (Figure 10), which is associated with the highest uncertainty measurements from the array of 48 GNSS stations, illustrates how the integrated dataset deviates significantly from the GNSS data in the up direction. Motions in the east-west and north-south directions are slightly more constrained, and motion in the up-down direction is significantly transformed after combining the DInSAR and GNSS datasets together. The integration algorithm weighs the DInSAR time series more than the GNSS data based on the errors of the individual datasets (see Supplemental Material).
Finally, the BLBP GNSS site ( Figure 11) is situated over an area of comparatively stable ground along the southwestern part of the island, where relatively low activity was maintained over the duration of our time series (green diamond, Figure 1). Leading up to the 2018 volcanic eruption, the ground at BLBP experienced approximately 1 cm of eastward motion, while remaining within ± 4 mm of motion in the north-south direction. At the time of the 2018 eruption, the BLBP station moved~12 cm to the east, remained steady in its north-south component of motion, and unclearly moved in the up-down direction. Nevertheless, over the entire time series, we detected that the ground at the BLBP station moved a total of~3 cm to the southeast, and~4 cm of subsidence occurred between November 2015 and April 2021.

Discussion
The 2018 event at Kilauea volcano was the largest caldera collapse and most effusive volcanic eruption in Hawaii within the past 200 years [26,[48][49][50][51]. Combining high-quality geodetic datasets, such as DInSAR and GNSS, resulted in detailed observations and precise measurements regarding the development and evolution of volcanic events.
The median uncertainty associated with the LOS DInSAR dataset was 6.60 mm, and the average was 6.54 mm. The maximum uncertainty associated with our LOS DInSAR displacement maps came from the final cumulative time slice and was equal to 1.03 cm. This uncertainty considers errors estimated from the DEM, orbital, and atmospheric corrections applied during processing.
The median errors associated with the raw GNSS data in the east, north, and up directions were 0.81 mm, 0.74 mm, and 3.27 mm, respectively. The corresponding maximum uncertainties associated with the raw GNSS data in the east, north, and up directions were 7.63 mm, 3.28 mm, and 9.65 mm. For direct comparison with the LOS DInSAR data, the converted LOS vector component of the raw GNSS data took on a maximum uncertainty of 1.27 cm.
Integrated uncertainties were combined from the raw DInSAR and GNSS sigma values using the same methodology as the real data (see Supplementary Material Equations S1-S8). Maximum integrated uncertainty in the LOS component of motion was 7.93 mm. When generating the 3D displacement maps, however, the error introduced from using the kriging algorithm and interpolating over undefined points in the spatial field must also be included, and in doing so, resulting in increased maximum uncertainties in the east, north, and up directions of 7.07 cm, 5.82 cm, and 5.71 cm, respectively.
At a single pixel, where kriging interpolation was not applied, the accuracy of our deformation measurements improved. For example, at the corresponding pixel where the CRIM GNSS station is situated, the median uncertainties for the raw GNSS time series data in the east, north, and up directions of motion were 0.77 mm, 0.70 mm, and 3.06 mm, respectively. The raw LOS DInSAR time series at this specific location had a median uncertainty of 1.99 mm. By integrating the two geodetic datasets together, the fused dataset had an improved median uncertainty of 0.76 mm, 0.70 mm, and 1.95 mm in the east, north, and up directions. Further analyses of the errors associated with each dataset are available in the Supplementary Material. Many other studies have utilized geodetic remote sensing techniques to better understand magmatic systems and volcanic eruptive history. Our raw interferometric products and LOS displacement values are consistent with those produced by the USGS HVO and European Space Agency, also processed with GMTSAR software with Snaphu phase unwrapping [52]. Furthermore, our results extend further evidence of the 62 summit collapse events between May 16 and August 2, reported by HVO. Another study by [50] used raw GNSS and DInSAR measurements from Hawaii as input for multi-reservoir conceptual models to quantitatively constrain the hydraulic connectivity between magmatic systems. By exploiting high-quality geodetic data, they learned that the Halema'uma'u magmatic reservoir is distinct from the South Caldera reservoir and that the ERZ is being fed simultaneously by both chambers. Our LOS inflation and deflation patterns over Kilauea and the ERZ agreed with their observations derived from descending interferograms from November 2018 through March 2019. Lundgren et al. [51] also used airborne InSAR measurements from the Glacier and Ice Surface Topography Interferometer (GLISTIN-A) instrument over Hawaii to measure the bulk volume of subaerial lava flows between 18 May and 15 September 2018. They found that 0.593 +/−0.011 km 3 sourced from the Lower ERZ and −0.836 +/−0.002 km 3 of material resulted from the summit collapse. Finally, a study at Piton de la Fournaise, France, used four interferograms to determine the displacement source and a temporal inversion of GNSS data to describe the dynamics of magmatic propagation [53].
These examples show the broad, interdisciplinary applications that come from a mix of DInSAR and GNSS monitoring over volcanic regions and further support the relevance and benefits of this study. A more accurate, fused dataset creates potential for improved volumetric analysis during large events, enhanced volcanic source modeling and mapping of magma chamber geometries or subsurface transport channels, and may provide a better means to forecast initial eruption sites or potential lava flow pathways along the surface.

Conclusions and Upcoming Work
Geodetic datasets such as DInSAR and GNSS time series contain valuable information for earth scientists. These data are particularly useful in analyzing the long-term evolution of volcanic systems. When integrated together, they deliver improved, more constrained estimates of 3D motions of the Earth's surface. Inflation estimates such as that detected in our integrated time series prior to both eruptions are useful for hazard response purposes and applications and underline the importance of deformation monitoring in volcanic regions. The methods performed here can be applied to other regions in the world and can be easily adapted to different volcanic systems.
The integrated 3D displacement maps and associated time series provided new estimates and details of the spatial and temporal dynamics of the 2018 Hawaiian eruption. In conjunction with seismic, tide gauge, gas emission, and thermal datasets, the individual InSAR and GNSS time series and the combined results presented here will be streamed through machine learning algorithms capable of identifying pixels or clusters of pixels that exhibit anomalous movement. This will facilitate the detection and analysis of precursor motion and can be used to inform early warning systems. Eventually, the system will run continuously, providing researchers with large streams of historic or near real-time data from an array of geodetic sensors.
With a continuous, scalable, near real-time data streaming architecture in place, the combined data will contribute to the early detection of major geophysical events such as volcanic eruptions and improve our understanding of the fundamental processes involved over the locations of interest. The fused dataset can be used to improve volcanic modeling systems, using the long-term surface response to map the subsurface magmatic plumbing system. The results presented here also have additional detailed insight into earthquake dynamics, and further research may lead to the distinction or isolation of certain geophysical phenomena within the signal. Finally, the application of machine learning algorithms to these datasets will improve our understanding of the connections between these volcanoes or how activity at one location may have implications for potential activity at another and be critical inputs for early warning applications. These results will provide researchers with better methods in real-time or near real-time to monitor volcanic regions and evaluate resulting surface deformations with millimeter precision, ultimately advancing scientific discovery and benefitting efforts to further protect human life and property.