Particle image velocimetry (PIV), a method based on image cross-correlation, is widely used for obtaining velocity fields from time-series of images of deforming objects. Rather than instantaneous velocities, we are interested in reconstructing cumulative deformation, and use PIV-derived incremental displacements for this purpose. Our focus is on analogue models of tectonic processes, which can accumulate large deformation. Importantly, PIV provides incremental displacements during analogue model evolution in a spatial reference (Eulerian) frame, without the need for explicit markers in a model. We integrate the displacements in a material reference (Lagrangian) frame, such that displacements can be integrated to track the spatial accumulative deformation field as a function of time. To describe cumulative, finite deformation, various strain tensors have been developed, and we discuss what strain measure best describes large shape changes, as standard infinitesimal strain tensors no longer apply for large deformation. PIV or comparable techniques have become a common method to determine strain in analogue models. However, the qualitative interpretation of observed strain has remained problematic for complex settings. Hence, PIV-derived displacements have not been fully exploited before, as methods to qualitatively characterize cumulative, large strain have been lacking. Notably, in tectonic settings, different types of deformation-extension, shortening, strike-slip-can be superimposed. We demonstrate that when shape changes are described in terms of Hencky strains, a logarithmic strain measure, finite deformation can be qualitatively described based on the relative magnitude of the two principal Hencky strains. Thereby, our method introduces a physically meaningful classification of large 2-D strains. We show that our strain type classification method allows for accurate mapping of tectonic structures in analogue models of lithospheric deformation, and complements visual inspection of fault geometries. Our method can easily discern complex strike-slip shear zones, thrust faults and extensional structures and its evolution in time. Our newly developed software to compute deformation is freely available and can be used to post-process incremental displacements from PIV or similar autocorrelation methods. @en

Satellite gravimetry and altimetry measurements record gravity and elevation changes, respectively, which are useful for determining mass and volume change of the Antarctic Ice Sheet. Common methods employ products from regional climate modeling and firn modeling to aid interpretation and to link volume changes to mass changes. Estimating deterministic parameters over defined time periods is a conventional way to evaluate those changes. To overcome limitations of deterministic analyses with respect to time-variable signals, we have developed a state-space model framework. Therein, we jointly evaluate four mass and volume data sets by coupling of temporal signal variations. We identify long-term signals of ice drainage basins that are observed by the satellite gravimetry mission GRACE and several satellite altimetry missions from April 2002 until August 2016. The degree to which we can separate long-term and short-term variations strongly depends on the similarity of the mass and volume change time series. For the drainage system of the Pine Island Glacier (West Antarctica), our results show noticeable variations of the long-term trend with an acceleration of the contribution of ice dynamics to the mass balance from −11 ± 8 to −58 ± 8 Gt a−1. Our results in Dronning Maud Land (East Antarctica) show a positive long-term contribution to the mass balance at almost a constant rate. The presented approach can fit time-variable changes without artificial selection of periods of interest. Furthermore, because we only enforce common long-term time variations between mass and volume data, differences in mean trend rates help to uncover model discrepancies.@en

A devastating tsunami struck Palu Bay in the wake of the 28 September 2018 Mw = 7.5 Palu earthquake (Sulawesi, Indonesia). With a predominantly strike-slip mechanism, the question remains whether this unexpected tsunami was generated by the earthquake itself, or rather by earthquake-induced landslides. In this study we examine the tsunami potential of the co-seismic deformation. To this end, we present a novel geodetic data set of Global Positioning System and multiple Synthetic Aperture Radar-derived displacement fields to estimate a 3D co-seismic surface deformation field. The data reveal a number of fault bends, conforming to our interpretation of the tectonic setting as a transtensional basin. Using a Bayesian framework, we provide robust finite fault solutions of the co-seismic slip distribution, incorporating several scenarios of tectonically feasible fault orientations below the bay. These finite fault scenarios involve large co-seismic uplift (>2 m) below the bay due to thrusting on a restraining fault bend that connects the offshore continuation of two parallel onshore fault segments. With the co-seismic displacement estimates as input we simulate a number of tsunami cases. For most locations for which video-derived tsunami waveforms are available our models provide a qualitative fit to leading wave arrival times and polarity. The modeled tsunamis explain most of the observed runup. We conclude that co-seismic deformation was the main driver behind the tsunami that followed the Palu earthquake. Our unique geodetic data set constrains vertical motions of the sea floor, and sheds new light on the tsunamigenesis of strike-slip faults in transtensional basins.@en

One of the primary observational data sets of sea level is represented by the tide gauge record. We propose a new method to estimate variability on decadal time scales from tide gauge data by using a state space formulation, which couples the direct observations to a predefined state space model by using a Kalman filter. The model consists of a time-varying trend and seasonal cycle, and variability induced by several physical processes, such as wind, atmospheric pressure changes and teleconnection patterns. This model has two advantages over the classical least-squares method that uses regression to explain variations due to known processes: a seasonal cycle with time-varying phase and amplitude can be estimated, and the trend is allowed to vary over time. This time-varying trend consists of a secular trend and low-frequency variability that is not explained by any other term in the model. As a test case, we have used tide gauge data from stations around the North Sea over the period 1980-2013. We compare a model that only estimates a trend with two models that also remove intra-annual variability: one by means of time series of wind stress and sea level pressure, and one by using a two-dimensional hydrodynamic model. The last two models explain a large part of the variability, which significantly improves the accuracy of the estimated time-varying trend. The best results are obtained with the hydrodynamic model. We find a consistent low-frequency sea level signal in the North Sea, which can be linked to a steric signal over the northeastern part of the Atlantic.@en

Large-scale strike-slip faults are associated with significant strain partitioning in releasing/restraining bends and often display map-view curvatures ending in horse-tail geometries. Such faults are commonly associated with indentation tectonics, where shortening in front of indenters is transferred laterally to transpression, strike-slip and the formation of transtensional/extensional basins. We investigate how these structurally distinct domains are kinematically linked by the means of a crustal-scale analogue modelling approach where a deformable crust is moved against a stable and rigid indenter. The modelling demonstrates that the geometry of the indenter is the major controlling parameter driving strain partitioning and deformation transfer from thrusting and transpression to strike-slip and transtension, whereas the rotation of the mobile plate controls the opening of triangular shaped transtensional basins. Flow of the ductile crust leads to the distribution of deformation over a wider area, facilitating strike-slip splaying into transtension/extension behind the indenter. Our results show a very good correlation with the Moesian indentation in the Carpatho-Balkanides system of South-Eastern Europe, where strain is partitioned around the dextral Cerna and Timok strike-slip faults and transferred to thrusting in the Balkanides part of the Moesian indenter and to transtension/extension in the neighbouring South Carpathians.@en

Large-scale strike-slip faults are associated with significant strain partitioning in releasing/restraining bends and often display map-view curvatures ending in horse-tail geometries. Such faults are commonly associated with indentation tectonics, where shortening in front of indenters is transferred laterally to transpression, strike-slip and the formation of transtensional/extensional basins. We investigate how these structurally distinct domains are kinematically linked by the means of a crustal-scale analogue modelling approach where a deformable crust is moved against a stable and rigid indenter. The modelling demonstrates that the geometry of the indenter is the major controlling parameter driving strain partitioning and deformation transfer from thrusting and transpression to strike-slip and transtension, whereas the rotation of the mobile plate controls the opening of triangular shaped transtensional basins. Flow of the ductile crust leads to the distribution of deformation over a wider area, facilitating strike-slip splaying into transtension/extension behind the indenter. Our results show a very good correlation with the Moesian indentation in the Carpatho-Balkanides system of South-Eastern Europe, where strain is partitioned around the dextral Cerna and Timok strike-slip faults and transferred to thrusting in the Balkanides part of the Moesian indenter and to transtension/extension in the neighbouring South Carpathians.@en

Most of the seismic moment release of the complex earthquake sequence beneath the South Sandwich Islands occurred on the central part of the SS megathrust. Significant aftershock activity indicates that the central and southern megathrust was subsequently activated, i.e., where young South America lithosphere is subducted. Seismic activity thus seems to have been restricted by the lateral termination in the south of the SS Trench. Relatively little energy release occurred on the northern part of the megathrust. It was hypothesized by Govers and Wortel (2005) that here the South America slab breaks away from the surface part of the plate at the active STEP. Geochemical observations and earthquake P-axes orientations do not seem to agree with the hypothesis and we investigate the cause. We show results of new physical analog lab models that aim to elucidate what controls the geometry of the lithospheric STEP Fault. We study lithospheric tearing in the process of STEP evolution, which is dynamically driven by the buoyancy of the subducting slab. In our experiments, the lithosphere as well as asthenosphere are viscoelastic media in a free subduction setup. A stress-dependent rheology plays a major role in localization of strain in tearing processes of lithosphere such as slab break-off. The results show that the highly curved northern plate boundary is a STEP Fault following from lithospheric tearing at a depth of ~100km. This is a modification of the original STEP model of Govers and Wortel (2005). This is consistent with available observations along the northern Sandwich plate boundary, and likely exists in other STEP regions. The region’s largest recorded event, the 1929 Mw 8.3 earthquake, may reflect horizontal extension perpendicular to the STEP fault, which is also expected based on our experiments.@en

In this study, we examine the effect of transient mantle creep on the prediction of glacial isostatic adjustment (GIA) signals. Specifically, we compare predictions of relative sea level (RSL) change from GIA from a set of Earth models in which transient creep parameters are varied in a simple Burgers model to a reference case with a Maxwell viscoelastic rheology. The model predictions are evaluated in two ways: first, relative to each other to quantify the effect of parameter variation, and second, for their ability to reproduce well-constrained sea level records from selected locations. Both the resolution and geographic location of the RSL observations determine whether the data can distinguish between model cases. Model predictions are most sensitive to the inclusion of transient mantle deformation in regions that are near-field and peripheral relative to former ice sheets. This sensitivity appears particularly true along the North American west coast in the region of the former Cordilleran Ice Sheet, which experienced rapid sea-level fall following deglaciation between 14 and 12 kyr BP. Relative to the Maxwell case, Burgers models better reproduce this rapid phase of regional postglacial sea-level fall. As well, computed goodness-of-fit values in this region show a clear preference for models where transient deformation is present in the whole or lower mantle, and for models where the rigidity of the Kelvin element is weakened relative to the rigidity of the Maxwell element. In contrast, model predictions of relative sea-level change in the far-field show weak sensitivity to the inclusion of transient deformation.@en

The Pannonian basin is an extensional back-arc basin that has undergone neotectonic inversion and is currently shortening. The understanding and quantification of present-day deformation processes during this inversion are still incomplete. To this end, we investigate the active deformation of the Circum-Pannonian region via the interpolation of GNSS-derived velocity field and the derivation of the strain rate fields. For the interpolation of the velocity field, we use ordinary kriging, a strochastic interpolation method. Our results show that estimating a strain rate field that is virtually free of short-wavelength noise requires the scaling of the velocity uncertainties, i.e. assuming a minimum standard deviation of 1 mm/yr in our case. The deformation of the Circum-Pannonian region is defined by the 2–3 mm/yr, NNE-directed motion of the Dinarides, and by the 0.5–1.5 mm/yr, WSW to SSW directed motion of the eastern areas (European foreland, East Carpathians, South Carpathians, Transylvanian basin). These opposite-sense motions define a large-scale, on average NE-SW shortening and transpression-type deformation in the Dinarides as well as in the Pannonian basin, while the East and South Carpathians undergo regional N–S extension. Neotectonic structures generally show good agreement with the strain rate field, for example in the Dinarides, Eastern Alps, or in the western Pannonian basin. However, the presence of fault-parallel shortening or biaxial shortening along sinistral neotectonic structures in the central and eastern Pannonian basin show some discrepancy between current geodetic and observed neotectonic deformation. The vertical velocity field shows dominantly 100 and 1000 km wavelength signals, the former is probably related to the response of the Pannonian lithosphere-asthenosphere system to neotectonic basin inversion, while latter can possibly be explained by far-field subsidence patterns induced by the mantle response to melting of the Fennoscandian ice sheet during the current interglacial period.@en

The Gravity Recovery and Climate Experiment (GRACE) mission (launched 2002) and the Gravity Field and Steady-State Ocean Circulation Explorer (GOCE) mission (March 2009 to November 2013) collected spaceborne gravity data for the preseismic and postseismic periods of the 2011 Tohoku-Oki earthquake. In addition, the dense Japan GeoNet Global Navigation Satellite Systems (GNSS) network measured with approximately 1050 stations the coseismic and postseismic surface displacements. We use a novel combination of GNSS, GRACE, and GOCE observations for a distributed fault slip model addressing the issues with gravimetric and geometric change over consistent time windows. Our model integrates the coseismic and postseismic effects as we include GOCE observations averaged over a 2 year interval, but their inclusion reveals the gravity change with unprecedented spatial accuracy. The gravity gradient grid, evaluated at GOCE orbit height of 265 km, has an estimated formal error of 0.20 mE which provides sensitivity to the mainly coseismic and integrated postseismic-induced gravity gradient signal of -1.03 mE. We show that an increased resolution of the gravity change provides valuable information, with GOCE gravity gradient observations sensitive to a more focused slip distribution in contrast to the filtered GRACE equivalent. The 2 year averaging window of the observations makes it important to incorporate estimates of the variance/covariance of unmodeled processes in the inversion. The GNSS and GRACE/GOCE combined model shows a slip pattern with 20 m peak slip at the trench. The total gravity change (≈200 μGal) and the spatial mapping accuracy would have been considerably lower by omitting the GOCE-derived fine-scale gravity field information.@en

Earthquakes involve mass redistribution within the solid Earth and the ocean, and as a result, perturb the Earth's gravitational field. For most of the shallow (<60 km) earthquakes with Mw > 8.0, the GRACE satellite gravity measurements suggest considerable volumetric disturbance of rocks. At a spatial scale of hundreds of km, the effect of volumetric change exceeds gravity change by vertical deformation; for example, negative gravity anomalies associated with volumetric expansion are characteristic patterns after shallow thrust events. In this study, however, we report contrasting observations of gravity change from two intermediate-depth (100–150 km) earthquakes of 2016 & 2017 Mw 8.0 (two combined) Papua New Guinea thrust faulting events and 2019 Mw 8.0 Peru normal faulting and highlight the importance of compressibility in earthquake deformation. The combined 2016/17 thrust events resulted in a positive gravity anomaly of 5–6 microGal around the epicenter, while the 2019 normal faulting produced a negative gravity anomaly of 3–4 microGal. Our modeling found that these gravity changes are manifestation of vertical deformation with limited volumetric change, distinct from gravity changes after the shallow earthquakes. The stronger resistance of rocks to volume change at intermediate-depth results in largely incompressible deformation and thus in a gravity change dominated by vertical deformation. In addition, malleable rocks under high pressure and temperature at depth facilitated substantial afterslip and/or fast viscoelastic relaxation causing additional vertical deformation and gravity change equivalent to the coseismic change. For the Papua New Guinea events, this means that postseismic relaxation enhanced coseismic uplift and relative sea level decrease.@en

In regions undergoing glacial isostatic adjustment present-day horizontal surface motion is observed to point mostly, but not always, away from the former ice load. To interpret these observations, we investigate the direction of horizontal velocities using glacial isostatic adjustment models. The direction is controlled by the opposing actions of inward mantle flow and outward lithosphere motion. In contrast with the prevailing idea that glacial isostatic adjustment-induced horizontal velocities point outward, we show that velocities can be either outward or inward. Immediately after deglaciation velocities point inward but change direction to outward after a time that is controlled by mantle viscosity. Present-day horizontal velocities point outward for a uniform mantle viscosity below 1020 Pa s and inward for above 1022 Pa s, with a combination of outward and inward in between. Our results help to interpret GPS-observed horizontal velocities in areas with varying mantle viscosity.@en

Stress-dependent nonlinear upper mantle rheology has a firm base in rock mechanical tests, where this nonlinearity results from dislocation creep of minerals. In the last few decades there has been some attention to nonlinear, power-law, materials for application in scaled analogue experiments for tectonic processes. However, studies describing the rheology of analogue materials with the same nonlinear dependency on stress as observed for lithospheric mantle materials at relevant stress levels, are still lacking. In this study we have developed and rheologically tested materials based on combinations of silicone polymers and plasticine, with the aim of obtaining a material that can serve as a laboratory analogue to the power-law rheology of olivine aggregates at lithospheric mantle conditions. From our steady-state creep tests we find that it is possible to obtain such a power-law material, with effective viscosities over relevant model stress ranges [5–4000 Pa] that allow for nonlinear deformation at laboratory time scales. We apply the developed material to a process where localized deformation of the lithosphere can be expected: slab break-off. We study this process using analogue models, where we apply the new nonlinear material to the lithospheric mantle domains, while we use Newtonian glucose to represent the low viscous asthenosphere. Now that we properly manage power-law behavior in our analogue lithosphere materials, we are able to model localized lithospheric tearing.@en

The 2011 Tohoku-Oki earthquake with 9.0 Mw led to an enormous mass redistribution originated from large deformation due to faulting and had a massive impact on the coastal area of eastern Japan. While the satellite gravity mission GRACE (Gravity Recovery and Climate Experiment) can detect the gravitational change caused by this tremendous event, slip distributions are usually derived from GPS, seismic and (in the more particular case) tsunami data. We evaluate the differences between measured and modeled coseismic gravity changes for three fault slip models derived from either GPS and tsunami data, GRACE data, or a combination of all three data types. The data are weighted according to their measurement accuracy in a Bayesian joint inversion approach. We perform a long term average of GRACE data, which increases sensitivity and reduces artefacts, and find that the postseismic gravity change leaks into the derived mean gravity field. We try to reduce this problem by averaging only 6 months of postseismic GRACE data, where the postseismic gravity signal, which superimposes onto the coseismic signal of ≈ 6μGal (for a geometric based model) peaks approximately 3 months after earthquake occurrence. Consequently fault slip models merely derived from GPS (10 days avg.) and tsunami data (<5h time span) show deviations of ≈ 2 μGal to a GRACE 6 monthly averaged combined solution which indicates the difference accumulated from the geometric and gravimetric modeling and the postseismic gravity signal in the GRACE data.@en

The 2011 Tohoku-Oki earthquake with 9.0 Mw led to an enormous mass redistribution originated from large deformation due to faulting and had a massive impact on the coastal area of eastern Japan. While the satellite gravity mission GRACE (Gravity Recovery and Climate Experiment) can detect the gravitational change caused by this tremendous event, slip distributions are usually derived from GPS, seismic and (in the more particular case) tsunami data. We evaluate the differences between measured and modeled coseismic gravity changes for three fault slip models derived from either GPS and tsunami data, GRACE data, or a combination of all three data types. The data are weighted according to their measurement accuracy in a Bayesian joint inversion approach. We perform a long term average of GRACE data, which increases sensitivity and reduces artefacts, and find that the postseismic gravity change leaks into the derived mean gravity field. We try to reduce this problem by averaging only 6 months of postseismic GRACE data, where the postseismic gravity signal, which superimposes onto the coseismic signal of ≈ 6μGal (for a geometric based model) peaks approximately 3 months after earthquake occurrence. Consequently fault slip models merely derived from GPS (10 days avg.) and tsunami data (<5h time span) show deviations of ≈ 2 μGal to a GRACE 6 monthly averaged combined solution which indicates the difference accumulated from the geometric and gravimetric modeling and the postseismic gravity signal in the GRACE data.@en

We aim to better understand the overriding plate deformation during the megathrust earthquake cycle. We estimate the spatial patterns of interseismic GNSS velocities in South America, Southeast Asia and northern Japan and the associated uncertainties due to variations in network density and observation uncertainties. Interseismic velocities with respect to the overriding plate generally decrease with distance from the trench with a steep gradient up to a ‘hurdle’, beyond which the gradient is distinctly lower and velocities are small. The hurdle is located 500–1000 km away from the trench for the trench-perpendicular velocity component, and either at the same distance or closer for the trench-parallel component. Significant coseismic displacements were observed beyond these hurdles during the 2010 Maule, 2004 Sumatra–Andaman, and 2011 Tohoku earthquakes. We hypothesize that both the interseismic hurdle and the coseismic response result from a mechanical contrast in the overriding plate. We test our hypothesis using physically consistent, generic, 3-D finite element models of the earthquake cycle. Our models show a response similar to the interseismic and coseismic observations for a compliant near-trench overriding plate and an at least five times stiffer overriding plate beyond the contrast. The model results suggest that hurdles are more prominently expressed in observations near strongly locked megathrusts. Previous studies inferred major tectonic or geological boundaries and seismological contrasts located close to the observed hurdles in the studied overriding plates. The compliance contrast probably results from thermal, compositional and thickness contrasts and might cause the observed focusing of smaller-scale deformation like backthrusting.@en

Before geodetically derived strain and rotation rates can be robustly compared to geological or seismological observations, we need reliable strain rate uncertainties. Various methods exist to compute strain rates from GNSS-derived interseismic velocities, but a realistic representation of interpolation uncertainties has remained a challenge. The main problem is that commonly used deterministic interpolation methods do not account for uncertainty resulting from the absence of information in between observation sites. We apply stochastic interpolation by means of ordinary kriging to propagate errors both from discontinuous data coverage as well as from observation uncertainties to our strain rate estimates. However, interseismic horizontal surface velocities in tectonically active regions are spatially highly non-stationary, with high spatial variability around active faults and lower velocity variability in tectonically more stable regions. This requires an extension of traditional ordinary kriging approaches. For interpolation uncertainties that reflect the local variability and spatial correlation of the observed surface velocities, we apply a novel method that incorporates the spatially variable statistics of the underlying data. We estimate realistic uncertainties and covariances of the interpolated velocity field. For regions with a high spatial velocity variability, we find a large increase in uncertainty with increasing distance from observation sites, while in areas with little spatial variability, we estimate a small increase in uncertainty with distance. Subsequently, we propagate interpolated velocity covariance to strain rate uncertainties, such that we can assess the statistical significance of the interpolated strain rate field. Applied to a number of actively deforming regions, including the Pacific coast of North America and Japan, we show to what degree we can robustly determine strain rates based on available GNSS-derived velocities. Realistic uncertainties assist the community to better discriminate continuous or localized deformation on active faults from the available geodetic data.@en

The GRACE and GRACE Follow-On era is a rare opportunity from which we can study the Earth’s response to great earthquakes across diverse tectonic settings at time scales from days to decades. Earthquakes cause abrupt gravity field changes and triggered postseismic adjustment expected to continue for years to decades. There were a total ~20 great earthquakes (Mw≥8.0)during the GRACE era (2002-2017). GRACE data has detected regional-scale coseismic and postseismic gravity changes, including 2004 Mw9.2 Sumatra-Andaman, 2005 Mw8.5 Nias, 2006/2007 Mw8.3 Kuril, 2007 Mw8.5 Bengkulu, 2009 Mw8.1 Samoa-Tonga, 2010 Mw8.8 Maule, 2011 Mw9.0 Tohoku-Oki, 2012 Mw8.6 Indian Ocean (Wharton Basin), and 2013 Mw8.3 Okhotsk earthquakes. These events include shallow megathrust, strike-slip, normal faults ruptures as well as a deep (~600 km) source earthquake. The GRACE gravity data demonstrated the importance of compressibility in governing large-scale coseismic gravimetric change. For the largest of these events (Mw≥8.6), we were able to constrain the seismic source depth and long-wavelength slip distribution, bulk modulus, elastic lithosphere thickness, sea water redistribution, and the Earth’s rheological structure. Most of the 20 events with Mw≥8.0, during the GRACE era, exhibited a postseismic response best approximated by bi-viscous relaxation with range of asthenospheric viscosities including a lower range of 5x10^17–10^18 Pa s. In this presentation, we review these events using monthly time series of the entire GRACE geopotential change datasets from 2002 to 2017. Using the elastic and viscoelastic normal mode computations, we analyze the GRACE data to optimally constrain the earthquake source models and the regional rheology models using a self-consistent approach. We determine the coseismic and postseismic geopotential change models by inverting the GRACE L2 data for those events. Additionally, our modeling results will be provided to the GRACE Project and science community so that they can be used to remove the earthquake-related GRACE response to enable climate investigations that rely on GRACE data free of tectonic signals. Our comprehensive products can be used to improve quantification of secular trends of ocean, cryosphere, and hydrological mass transport from GRACE and GRACE Follow-On measurements. @en

The GRACE and GRACE Follow-On era is a rare opportunity from which we can study the Earth’s response to great earthquakes across diverse tectonic settings at time scales from days to decades. Earthquakes cause abrupt gravity field changes and triggered postseismic adjustment expected to continue for years to decades. There were a total ~20 great earthquakes (Mw≥8.0)during the GRACE era (2002-2017). GRACE data has detected regional-scale coseismic and postseismic gravity changes, including 2004 Mw9.2 Sumatra-Andaman, 2005 Mw8.5 Nias, 2006/2007 Mw8.3 Kuril, 2007 Mw8.5 Bengkulu, 2009 Mw8.1 Samoa-Tonga, 2010 Mw8.8 Maule, 2011 Mw9.0 Tohoku-Oki, 2012 Mw8.6 Indian Ocean (Wharton Basin), and 2013 Mw8.3 Okhotsk earthquakes. These events include shallow megathrust, strike-slip, normal faults ruptures as well as a deep (~600 km) source earthquake. The GRACE gravity data demonstrated the importance of compressibility in governing large-scale coseismic gravimetric change. For the largest of these events (Mw≥8.6), we were able to constrain the seismic source depth and long-wavelength slip distribution, bulk modulus, elastic lithosphere thickness, sea water redistribution, and the Earth’s rheological structure. Most of the 20 events with Mw≥8.0, during the GRACE era, exhibited a postseismic response best approximated by bi-viscous relaxation with range of asthenospheric viscosities including a lower range of 5x10^17–10^18 Pa s. In this presentation, we review these events using monthly time series of the entire GRACE geopotential change datasets from 2002 to 2017. Using the elastic and viscoelastic normal mode computations, we analyze the GRACE data to optimally constrain the earthquake source models and the regional rheology models using a self-consistent approach. We determine the coseismic and postseismic geopotential change models by inverting the GRACE L2 data for those events. Additionally, our modeling results will be provided to the GRACE Project and science community so that they can be used to remove the earthquake-related GRACE response to enable climate investigations that rely on GRACE data free of tectonic signals. Our comprehensive products can be used to improve quantification of secular trends of ocean, cryosphere, and hydrological mass transport from GRACE and GRACE Follow-On measurements. @en