Analysis of GRACE Satellite Data in Terms of Geomonitoring of Strong Underwater Earthquakes Valentin Kashkin2, Tatyana Rubleva1 and Konstantin Simonov2 1 Siberian Federal University, 79 Svobodny st., Krasnoyarsk 660041, Russia 2 University 1, Address, City, Index, Country Institute of Computational Modelling of the Siberian Branch of the Russian Academy of Sciences, 50/44 Akademgorodok, 660036, Krasnoyarsk, Russia Abstract The study is devoted to a detailed analysis of observational data from the GRACE space system to refine the parameters of the focal zone of the 2010 catastrophic earthquake with the magnitude MW = 8.8 which occurred in the subduction zone near the Chilean coast. A technique was developed for con-structing maps of the distribution of the EWH parameter (equivalent water height above the geoid contour) and interpreting their anomalies relative to the source area. It is shown that there is a negative correlation between the value of the geodynamic parameter H (the distance from the hypocenter to the position of the Earth-Moon barycenter) for this strong earthquake with МW = 8.8 and the parameter EWH. Within the framework of the comparative analysis, digital maps of changes in the EWH parameter were also constructed for the studied focal areas of the Chilean earthquakes of 2014 and 2015. Keywords 1 Gravitational field, EWH parameter, space systems GRACE, satellite data processing, barycenter of the moon-earth system, correlation analysis 1. Introduction The temporal variations of the local gravity field observed during strong earthquakes were studied [1, 2] on the basis of gravimetric ground-based and satellite data significantly supplemented by the new data on the study of the Earth's gravitational field and its anomalies owing to the projects implemented in the 21st century by NASA-DLR space agencies, such as CHAMP (Challenging Mini- satellite Payload), GOCE (Gravity Field and Steady-State Ocean Circulation Explorer), GRACE (Gravity Recovery and Climate Experiment) and GRACE-FO (Gravity Recovery and Climate Experiment Follow-On). In [1] the duration of the period of gravitational changes during the seismic activity in Indonesia in December 2004 (earthquake with the magnitude Mw = 9.0) was revealed. The period of such variations was approximately 28 days. The time interval when significant gravitational changes were observed was 23 days. This is the time of activation of the foreshock process in this seismically active area. In [2] based on the data of the space system (CS) GRACE, estimates were made of the possibility of identifying coseismic variations in the local gravitational field during strong earthquakes with a magnitude MW of more than 8.0. It was determined that the probability of detecting signals from earthquakes with the magnitude Mw = 9.0 was higher than 98%, with Mw = 8.8 it was about 60%, and with Mw = 8.6 about 40%. In [3] it is hypothesized that variations in the gravitational field reflect the hydrological movements of water masses around the solid Earth which are caused by an increase in the elastic surface load on faults. The EWH (Equivalent Water Height) parameter or the equivalent water height SibDATA 2021: The 2nd Siberian Scientific Workshop on Data Analysis Technologies with Applications 2021, June 25, 2021, Krasnoyarsk, Russia EMAIL: simonovkv@icm.krasn.ru (K. Simonov); rtcvbk@rambler.ru (V. Kashkin); Tvrubleva@mail.ru (T. Rubleva) ORCID: 0000-0001-2345-6789 (V. Kashkin); 0000-0003-1876-2858 (T. Rubleva); 0000-0002-6829-3087 (K. Simonov) © 2021 Copyright for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0). CEUR Workshop Proceedings (CEUR-WS.org) above the geoid contour (in cm) [4] is used as an amplitude of the gravity change in the GRACE and GRACE-FO data. The geoid in [5] is understood as an equipotential surface (balanced) surface of the Earth's gravitational field which approximately corresponds to the topographic features of the planet and level of the World Ocean in an undisturbed state. The purpose of our study is to search for patterns in the distribution of the EWH parameter in the region of focal zones of the strongest underwater earthquakes during the activation period of the foreshock process. The retrospective analysis of geophysical information about a series of three strongest earthquakes with the magnitude МW> 8.0 which occurred in a single subduction zone along the Pacific coast of Chile in 2010, 2014, and 2015 was performed. To clarify the geodynamic and gravitational parameters during the preparation of these catastrophic underwater earthquakes, a method for comparative analysis of the GRACE measurements and data from the USGS global geodynamic monitoring database was developed and tested. The studies show that the identification of the features in the gravitational field of the Chilean region and the comparative analysis of the main geophysical parameters obtained using the measurement information of two observation systems (GRACE CS and USGS geomonitoring) in the studied source zones allows solving urgent prognostic problems. 2. Application of the GRACE technology for solving geomonitoring problems At present, a large amount of data has been accumulated from the space systems GRACE (2002- 2017) and GRACE-FO (from 2018 to the present) [6]. On the basis of this information, the space-time variations of the Earth's gravity field with a period of 30 days are studied [7, 8]. The GRACE CS implements the Satellite-to-Satellite Tracking (SST) method [8]. High-precision positioning of Grace 1 and Grace 2 is provided by linking to the GPS navigation system. The discrepancy between the val- ues of the binding with the ground data of the laser ranging is 2-3 cm [6]. The range and rate of change in the distance between the Grace 1 and Grace 2 satellites are determined at altitudes of 300- 500 km. The measurement error of the GRACE system does not exceed 10 microns (1 micron = 10-6 m). Note that the accuracy of the GRACE data was confirmed by numerous scientific studies [2, 9, 10]. Although this system measures the characteristics of the gravitational field over the water surface and land with the same accuracy in geodynamic problems, it is necessary to solve the problem of separating deep and near-surface effects. For this, it is proposed to select geophysical objects during the period of activation of geodynamic processes in the transition zones Land - Sea. In our study, we analyze the distribution of the EWH parameter in the source zones of strong earthquakes in the Chilean region. Its values according to [4] are determined taking into account satellite measurements and the calculated harmonic coefficients of the geopotential relative to the mean field of the geoid (model). The information on this parameter is given on the website [6]. On this basis, an archive of the GRACE satellite data was formed for the studied focal zones. 3. Computational technique for constructing and analyzing digital maps of the EWH parameter distribution In our research, to solve the set tasks of the retrospective analysis and interpretation of GRACE data in seismically active regions of strong underwater earthquakes, a computational technique was proposed including the following stages: analysis of variations in the EWH parameter in the disturbed geoenvironment relative to the epicenter area; construction of digital maps of the EWH parameter distribution under the disturbed (EWH2) and background (EWH1) conditions relative to the focal zone; construction of a digital map of the distribution of the relative deviations EWH determined from the difference EWH = EWH2 - EWH1 in the study area; detection and description of the EWH characteristic relative to the location of the focal zone and hypocenter of the studied strongest earthquake; calculation of the geodynamic parameter H defined as the distance from the hypocenter to the position of the barycenter of the Earth-Moon system; estimation of the correlation coefficients between the series of the H and EWH values; refinement of the spatio-temporal characteristics of the focal zones of strong underwater earthquakes with respect to the identified anomalies in the distribution of the EWH parameter. In [11] we studied the processes of preparation of the strongest earthquakes in seismically active regions of the Pacific Ocean based on the created computational technique and analysis of digital EWH maps. It was revealed that during the preparation and relaxation of the earthquakes the equivalent water height above the geoid in seismically active regions changed significantly. However, there remains the questions whether the EWH parameter is suitable for specifying the geodynamic features during the preparation of the strongest earthquake and whether this parameter can be used as a characteristic for identifying the relationship between the stress-strain state of the geomedium and foreshock processes. 4. Experimental studies of gravitational anomalies for the area of the strongest Chilean earthquakes The focal zone of the strongest tsunamigenic Chilean earthquake was studied where 217 post- seismic events with the magnitudes МW> 5 and hypocenter depths from 11 to 50 km were recorded during the year from February 27, 2010 to February 28, 2011 [12]. These aftershocks arose after the catastrophic earthquake with MW = 8.8 which occurred off the coast of Chile at 06:34:11 UTC and generated a strong tsunami. The focus of this seismic event was at a depth of 22.9 ± 9.2 km, its spatial dimensions were 600x120 km. The epicenter of the MW = 8.8 earthquake with the coordinates 36.122°S and 72.898°W was located in the ocean in the coastal zone of Central Chile (Maule region). According to the data on the geodynamics of the study area [13], the Nazca plate moves under the South American plate at a speed of 66 mm / year. The seismotectonic history of this region is described in [14]. The features of the for-mation of the source of the Chilean earthquake MW = 8.8 and the seismic wave field were analyzed based on the key model of the subduction zone in [15]. Based on the GRACE data a graph of variations in the EWH parameter in the disturbed geoenvironment relative to the epicenter area of the catastrophic earthquake MW= 8.8 was constructed for 25 months including the month of the main seismic event (Figure 1). Figure 1: Variations in the equivalent water level (EWH) relative to the epicenter area of the Maule earthquake МW = 8.8 This period covers 2 time intervals of the geodynamic activity: the first is the activation phase of the foreshock process (from January 2009 to January 2010), and the second is the phase of post- seismic events (from March 2010 to February 2011). The values "0" on the vertical scale of EWH (Fig. 1) correspond to the level when the surfaces of the geoid and waters of the World Ocean almost coincide. The variations of this parameter relative to the epicenter of the 2010 Chilean earthquake were estimated in the time intervals before and after the main moment of the earthquake with MW=8.8 in order to identify a "precursor" (Fig. 1). Note that negative EWH values were observed during the 7- month time interval from March to October 2009. The only exception was a slight increase in the EWH parameters in May and September. In May 2009, EWH = 0.2 cm and in September EWH was 0.44 cm. The minimum value of this parameter in February 2010 was equal to EWH = -2.58 cm. It should also be noted that during the recording of the postseismic events in the period from April 2010 to July 2011 in the studied region the EWH values were positive. Earlier in [16], according to the GRACE data, the monitoring of gravity anomalies in the studied region was carried out. Anomalous areas were observed after the main seismic event with МW = 8.8. A month after the Maule earthquake to the east of its epicenter, a negative gravity anomaly (–5 μGal) was recorded. Coseismic effects also revealed a smaller positive gravity anomaly. The analysis of the spatial structure of the observed anomalies confirmed the previously advanced theoretical assumptions on the redistribution of the internal mass of the solid Earth during the geodynamic activity. Based on the GRACE satellite data [6], we constructed digital maps of the EWH parameter distribution under the disturbed and background conditions relative to the source zone. The relative deviations EWH were determined from the difference between the parameters, EWH2 - EWH1 = EWH. Here, the EWH2 values show the distribution of the equivalent water level above the geoid in February 2010 in the disturbed geoenvironment. The EWH1 values are related to the source zone in the background geophysical environment. The analysis of the USGS geomonitoring catalog [12] shows that in February 2008 no seismic events were recorded in the future source of the earthquake with МW = 8.8. Taking these results into account, the spatial distribution of EWH1 was obtained under relatively quiet “background” conditions of the geoenvironment. It was found that in the epicenter area of the future Maule earthquake in February 2008 the value of EWH1 was 0.89 cm. The visualization of the distributions of the equivalent water level above the geoid relative to the location of the focal zone of the studied strongest earthquake is shown in Figure 2. Figure 2: Digital maps of the spatial values EWH2 (a), EWH1 (b) and deviations EWH (c) relative to the source zone of the earthquake with МW = 8.8 Here, the "cross" denotes the epicenter of the earthquake. Thus, the spatial location of the EWH2 values in the disturbed geoenvironment for February 2010 is shown in Figure 2 a. The image of the EWH1 background values under relatively quiet seismic conditions in February 2008 is given in Figure. 2 b. The type of the deviations EWH relative to the focal zone of the earthquake МW = 8.8 is shown in Figure 2 c. In Figure 2 a and Figure 2 c it can be seen that the anomalous region with the negative EWH values (negative anomaly) is closer to the equatorial latitudes. It is oriented towards SW-NE. The epicenter of the February 27, the 2010 earthquake is located on the periphery of this anomaly in the middle latitudes of the Southern Hemisphere. In Figure 2, a positive anomalous region is observed to the west of the negative anomaly and to the south of the earthquake epicenter. Thus, our results confirmed the results obtained earlier, thus evidencing the redistribution of the internal mass of the solid Earth during the strongest earthquakes. In addition, the constructed differential digital map of the distribution of the deviations EWH (Fig. 2 c) made it possible to clarify the location of the earthquake hypocenter relative to two spatial regions (negative and positive anomalies in terms of the EWH values). To study the influence of external geophysical processes, the geodynamic parameter H was calculated. The parameter H is the distance from the hypocenter of the 2010 catastrophic earthquake (MW = 8.8) to the Earth-Moon barycenter in the period from January 2009 to December 2011 (Fig. 3). Figure 3: Variations in the EWH and H parameters for the 2010 Chilean earthquake The method for calculating the coordinates of the barycenter is described in [17]. Correlation analysis revealed a relationship between the geodynamic parameter H and the EWH characteristic. The largest correlation coefficient between the H and EWH series had a value of R = -0.61 in the period from January 2010 to January 2011. At this time the maximum number of postseismic events is observed in the Central Chile region. The next studied strong earthquake with МW = 8.2 occurred on April 1, 2014 at 23:46:47 UTC in the coastal zone of Central Chile (in the region of the Iquique province) and generated a noticeable tsunami. The coordinates of the epicenter were 19.610°S and 70.769°W. The outbreak was located at a depth of 25 km. In the period from April 1 to April 15, 2014 according to the USGS data [12], there occurred 376 aftershocks. According to the developed methodology, the values of the EWH parameter were estimated relative to the hypocenter of this earthquake before and after the main moment of the earthquake with МW = 8.2 from April 2013 to April 2015 in order to identify the "precursor". The analysis of the variations in the equivalent water height above the geoid made it possible to reveal a 7-month time interval of the negative EWH values during the preparation of the studied seismic event. This was observed in June - December 2013. Noteworthy is the fact that during the activation of the foreshock process, January - February 2014, the EWH values were positive. The minimum value of this parameter in April 2014 was EWH = -1.1 cm. It should also be noted that during the recording of the postseismic events from May to June 2014 there was a significant increase in EWH in the epicenter area of the studied region. This analysis is complemented by the construction of a digital map of the distribution of the EWH1 parameters in the disturbed geoenvironment and in the background conditions. Based on the geodynamic information from the USGS catalog [12], April 2008 was chosen as the background geoenvironment. This time can be considered the most seismically calm for the region under study. In the disturbed conditions of the geoenvironment in April 2014 two negative and one positive anomalies were observed. The hypocenter of the strongest earthquake in the Iquique region was located north of the center of the first negative anomaly in equatorial latitudes. An examination of the spatial distribution of the EWH1 parameter revealed the orientation of this anomaly. It is directed from the southwest to the northeast. Based on the GRACE satellite data [6], a difference digital map of the distribution of the parameter EWH1 relative to the focal zone of this earthquake with МW = 8.2 was constructed. As in the previously considered case, the epicenter of this seismic event was located between two anomalous areas with the opposite EWH1 values. To the left of it, there is an uplift of the ocean floor (positive values of EWH1), and to the right there is an area of subsidence of the surface (negative values of EWH1). Consequently, the observed temporal changes in the gravity field indicate a redistribution of the internal mass in the lithosphere. An earthquake with МW = 8.3 occurred on September 16, 2015 at 22:54:32 UTC in the coastal zone of Central Chile (in the region of Ilhapel province) and as the previous ones it formed a significant tsunami. The coordinates of its epicenter were 31.573°S and 71.674°W, and the source depth was 22.4 km [12]. For the peak region of this seismic event according to the GRACE measurements, a detailed analysis of the distribution of the EWH2 values in the disturbed geoenvironment was made. Taking into account the methodology, the values of the EWH parameter relative to the hypocenter were estimated before and after the main moment of the earthquake with MW = 8.3. It was revealed that all the observed variations of this parameter during the study period September 2014 - June 2016 were positive. The minimum value of EWH2 at the epicenter (September 2015) was 0.78 cm. Note that during the activation of foreshock processes, April - August of this year, the EWH2 values were significantly higher than the minimum. From October 2015 to June 2016 inclusive, an intense increase in postseismic effects was recorded relative to the source of the Chilean earthquake. For the investigated seismically active region, we built a set of digital maps for the disturbed and background geoenvironment as well as an image of the relative deviations relative to the peak area of the earthquake MW = 8.3. It was revealed that the hypocenter of the earthquake with MW = 8.3 was located in the western part of the periphery of the positive anomaly. This anomalous area was oriented southeast (SE) - northwest (NW). To the left of the hypocenter there is an anomaly with the positive EWH2 values. This is where the ocean floor rises. And to the right there is an area with the negative values of EWH2 where the subsidence of the surface takes place. To interpret the obtained results, the following series of observations were made. The EWH series characterizes the variations in the EWH parameter at the epicenter of the catastrophic Chilean earthquake of 2010 with МW = 8.8 (from August 2009 to August 2010), the EWH1 series were obtained for the 2014 earthquake with МW = 8.2 (from October 2013 to October 2014) as well as the EWH2 series for the 2015 earthquake with МW = 8.3 (from March 2015 to March 2016). The correlation coefficients between these time series were estimated. The highest correlation coefficients equal to 0.68, 0.78, and 0.98 between the above series were to be characteristic for a period of 4 months, including 3 months of preparation for the studied earthquake and 1 month when the seismic event itself occurred. Thus, using the correlation analysis, the temporal variations of EWH in the series of the Chilean seismic events were found to be directly related to the preparation of strong earthquakes with a magnitude of M higher than 8.0. 5. Conclusion Based on the developed computational methodology for the analysis and interpretation of the GRACE measurement data and global seismic monitoring data the focal area of the catastrophic Chilean earthquake of 2010 with the magnitude MW = 8.8 was studied in comparison with the subsequent strongest earthquakes in this region which occurred in 2014 and 2015. Digital maps of changes in the EWH parameter over the geoid contour were constructed for the investigated seismically active area of these Chilean earthquakes as well as for the “background states” of the seismic regime of this area. A negative anomaly in the EWH parameter was identified which was observed in February 2010 to the west relative to the hypocenter of the 2010 catastrophic earthquake and oriented in the NE-SW direction. In 2014 and 2015, anomalous regions in the EWH parameter were also observed in the seismically active regions of these strong earthquakes. We believe that the formation of the corresponding anomalous area in terms of the EWH parameter is associated with geodynamic processes during this observation period. It is also shown based on the calculation results that there is a negative correlation between the series of the geodynamic parameters H and EWH for the catastrophic Chilean earthquake of 2010 with МW = 8.8. The highest correlation coefficient between the H and EWH series, equal to R = 0.61 was obtained when recording the maximum number of the postseismic events in the studied region. Using the correlation analysis, the temporal variations of EWH in the series of Chilean seismic events was found to be directly related to the preparation of strong earthquakes with the magnitude of M higher than 8.0. In the course of further research, it is planned to study a representative class of the strongest earthquakes in the Pacific region for a detailed analysis of the peculiarities of the variability of the EWH parameter and its relationship with global geodynamic and geophysical processes. 6. References [1] V. E. Khain, E. N. Khalilov, Gravitational effects before strong remote earthquakes. Bulletin of the International Academy of Sciences 2 (2007) 46–52. [2] V. O. Mikhailova, I. Panetc, M. Hayn, E. P. Timoshkina, S. Bonvalot, V. Lyakhovsky and et al.: Comparative study of temporal variations in the earth’s gravity field using GRACE gravity models in the regions of three recent giant earthquakes, Physics of the Earth. 5 (2014)171–191. [3] Y. Mitsui, K. Yamada, Possible correlation between annual gravity change and shallow background seismicity rate at subduction zone by surface load, Earth, Planets and Space 69 (2017) 1–7. [4] J. Wahr, M. Molenaar, F. Bryan,Time variability of the Earth's gravity field: Hydrological and oceanic effects and their possible detection using GRACE, J. Geoph. Res. Solid Earth. 103 (1998) 30205–30229. [5] B. Hoffmann-Wellengoff, G. Moritz, Physical Geodesy, Publishing house MIIGAiK, M, 2007. [6] NASA. GRACE, URL: https://grace.jpl.nasa.gov/. [7] B. D. Tapley, S. Bettadpur, J. C. Ries, P. F. Thompson, M. M. Watkins, GRACE measurements of mass variability in the Earth system, Science 305 (2004) 503–505. [8] A. Peidou, S. Pagiatakis, Gravity Gradiometry with GRACE Space Missions: New Opportunities for the Geosciences, Journal of Geophysical Research (Solid Earth) 124 (2019) 9130–9147. [9] V. O. Mikhailov, M. Diaman, A. A. Lyubushin, E. P. Timoshkina, S. A. Khairetdinov, Large- scale seismic creep in areas of strong earthquakes according to data from GRACE satellites on temporal variations of the gravitational field, Physics of the Earth. 5 (2016) 70–81. [10] N. S. Tkachenko, I. V. Lygin, Application of the GRACE satellite mission for solving geological and geographical problems, Bulletin of Moscow University, Series 4: Geology 2 (2017) 3–7. [11] K.V. Simonov, V. B. Kashkin, T. V. Rubleva, K. V. Krasnoshekov, Analysis of GRACE satellite measurements over seismically active areas of the strongest earthquakes, in: E3S Web of Conferences: Regional Problems of Earth Remote Sensing (RPERS 2018), volume 75, 2019, pp.1–5. [12] USGS, URL: https://earthquake.usgs.gov/. [13] M. Moreno, M. Rosenau, O. Oncken, Maule earthquake slip correlates with pre-seismic locking of Andean subduction zone, Nature 467 (2010) 198–202. [14] I. S. Vladimirova, Modeling of postseismic processes in subduction zones, Geodynamics and Tectonophysics 3 (2012) 167–178. [15] R. Kh. Mazova, Kh. F. Ramirez, N. A. Baranova, A. G. Rassadin, Catastrophic earthquakes and tsunamis in Chile. Evidence of a justified forecast, Proceedings of NSTU 2 (2014) 43–52. [16] S‐C. Han, J. Sauber, R. Riva, Contribution of satellite gravimetry to understanding seismic source processes of the 2011 Tohoku-Oki earthquake, Geophys. Res. Lett. 38 (2011)L24312. [17] V. G. Sibgatulin, S. A. Peretokin, A. A. Kabanov, Resonances of gravitational tides are a powerful energy source of geodynamic processes in the earth's crust, Journal of the Siberian Federal University. Technics and Technology. 9 (2016) 146–165.