Treatment of discontinuities inside Earth models: Effects on computed coseismic deformations

In this paper, we study how coseismic deformations calculated in 1066 Earth models are affected by how the models treat Earth discontinuities. From the results of applying models 1066A (continuous) and 1066B (discontinuous), we find that the difference in Love numbers of strike‐slip and horizontal tensile sources are bigger than dip‐slip and vertical tensile sources. Taken collectively, discontinuities have major effects on Green’s functions of four independent sources. For the near‐field coseismic deformations of the 2013 Okhotsk earthquake (Mw 8.3), the overall differences between theoretical calculations in vertical displacement, geoid, and gravity changes caused by discontinuities are 10.52 percent, 9.07 percent and 6.19 percent, with RMS errors of 0.624 mm, 0.029 mm, and 0.063 μGal, respectively. The difference in far‐field displacements is small, compared with GPS data, and we can neglect this effect. For the shallow earthquake, 2011 Tohoku‐Oki earthquake (Mw 9.0), the differences in near‐field displacements are 0.030 m (N‐S), 0.093 m (E‐W), and 0.025 m (up‐down) in our study area with the ARIA slip model, which gives results closer to GPS data than those from the USGS model. The difference in vertical displacements and gravity changes on the Earth’s surface caused by discontinuities are larger than 10 percent. The difference in the theoretical gravity changes at spatially fixed points truncated to degrees 60, as required by GRACE data, is 0.0016 μGal and the discrepancy is 11 percent, with the theoretical spatial gravity changes from 1066B closer to observations than from 1066A. The results show that an Earth model with discontinuities in the medium has a large effect on the calculated coseismic deformations.


Introduction
Seismic dislocation theory is an important basic theory in seismology and geodesy. It plays a key role in seismic fault inversion and geodetic observation data interpretation. Seismic dislocation theory is based on a regular geometric Earth model, such as a semiinfinite space, homogeneous sphere, or layered sphere. Okada (1985) gave a set of simple analytical formulas, based on previous research results, that are widely used to calculate the coseismic displacement and strain changes in the homogeneous elastic halfspace Earth model. Later, Okada (1992) extended this formulation to the calculation of seismic deformation in the Earth's interior. Okubo (1992) gave analytical formulas for calculating the gravita-tional potential and gravitational change of the medium model in a semi-infinite space. Because of its mathematical simplicity, dislocation theory based on semi-infinite space Earth models is still widely used. However, later researchers (Wang et al., 2003(Wang et al., , 2006 have improved the dislocation theory of the half-space model by taking into consideration the radial structure of the Earth, providing an effective scheme for calculating coseismic and postseismic deformation caused by any point source or generalized finite fault source. Sun and Okubo (1993) and Pollitz (1996) proposed a dislocation theory for a spherical elastic Earth model, which could calculate the seismic deformation over all of the Earth's surface. For the post-seismic relaxation computation, Pollitz et al. (2006) also considered a spherically-stratified, compressible, and self-gravitating viscoelastic Earth. Tanaka et al. (2006Tanaka et al. ( , 2007 introduced a new method to compute global post-seismic deformation in a spherically symmetric, self-gravitating, viscoelastic Earth model and investigated the effects of layering of the viscoelastic structure. At the same time, many researchers (Piersanti et al., , 1997Sabadini et al., 1995;Sabadini and Vermeersen, 1997;Soldati et al., 1998) introduced the normal model method to the incompressible spherical model with stratified structure, which further developed and improved seismic dislocation theory for the spherical Earth model.
Further efforts have then been carried out in the study of the influence of different Earth models on the calculation of seismic deformation in order to prove the importance of layered Earth models in practical applications. Pollitz (1996) investigated the coseismic deformation of radial layered structures and found that the influence of layered structures can reach 20%. Sun and Okubo (2002) found that the stratification and curvature effects on coseismic surface deformation (vertical displacement) can be larger than 20%. Dong J et al. (2014) found that the stratification effect for a point source simulation can be as high as 25%.
μ Basically, these studies were based on comparisons between homogeneous and inhomogeneous earth models. The radial structure of the Earth model is thus important for accurately calculating seismic deformation. On the slip model inversion, Hearn and Bürgmann (2005) studied the effect of elastic layering. They found that the recovered seismic moment is up to 40% greater for models incorporating depth-dependent shear modulus ( ) than it is for uniform elastic half-space models. However, these studies consider only the stratified structure of the Earth model, that is, the influence of the change in the average structure of the Earth model in the radial direction relative to that of a homogeneous sphere on the calculation of seismic deformation, and do not consider radially discontinuous interfaces (i.e., discontinuities in the medium parameters) in the Earth's interior. In the algorithms implementing their theoretical formulation, they discretize the radial layers of Earth model into a fine layering; the parameters around the discontinuities are then smoothed.
Within the philosophy that inspired the above listed studies, we herein make a step forward by considering the potential effects of Earth's discontinuities by comparing the results of seismological models against real data, from gravity and from surface displacements, from shallow and deep earthquakes, in the near and in the far fields. We thus report how the continuous or discontinuous structure of the Earth's interior medium affects the calculation of seismic deformation: this subject is in fact still an open question. Gilbert and Dziewonski (1975) proposed two Earth models, namely, 1066A (continuous model) and 1066B (discontinuous model), when they studied the structure of the Earth with the 1066 Earth model. Both of these models fit the original data very well. In the outer core and lower mantle, below a depth of 950 km, there is little difference between the two models. In practical applications, researchers usually use the 1066A model, rather than the 1066B model, without a self-consistent explanation for this preferential choice. Although both the 1066A and 1066B Earth models reflect the radially inhomogeneous structure of the Earth, how the difference between continuous and discontinuous structures affects the calculation of seismic deformation is not clear. To answer this question, the best approach is to examine the theoretical results within a dislocation theory based on different Earth models and give a statistical analysis of the discontinuity effects when compared with empirical data.
Therefore, in view of the structural differences between continuity and discontinuity in the 1066A and 1066B Earth models, this study estimates their effects on Love numbers and Green's functions (Sun WK et al., 2009) for the coseismic deformations caused by point sources. For a deep case study, we calculate the effect of coseismic deformation caused by the 2013 Okhotsk earthquake (M w 8.3); for a shallow case, we choose the 2011 Tohoku-Oki Earthquake (M w 9.0). We then discuss the effect of discontinuities, present empirical data and model predictions, and perform a statistical analysis of their similarities and differences.
2. Structure of 1066A and 1066B Earth Models μ As mentioned above, Gilbert and Dziewonski (1975) proposed two Earth models, namely, 1066A (continuous model) and 1066B (discontinuous model), whose medium model structures (Shear modulus ) are shown in Figure 1. The first discontinuous interface is defined at 11.0 km and 21.0 km, respectively. The two models are almost identical in the fluid outer core and in the lower mantle. There is little difference within the core, but differences are great within the upper mantle. The 1066A model (dashed line in Figure 1) assumes a smooth and continuous upper mantle,  Figure 1) features a discontinuous upper mantle with depths of discontinuity at 420 km and 671 km. Gilbert and Dziewonski (1975) assumed that the magnitudes of discontinuities are less than 0.08% and that an average distance of 200 km seems sufficient to eliminate the main differences between models 1066A and 1066B.
We thus evaluate the effects of these discontinuities on coseismic deformation caused by the point source and finite fault model, from the perspective of both surface displacements and gravity.

Effects of Discontinuities on Love Numbers and
Green's Functions for Point Sources Sun WK et al. (2009) proposed a set of numerical integration formulas to calculate coseismic deformation based on the spheric-h n , l n , k n , l t n ally symmetric, non-rotating, perfectly elastic, isotropic (SNREI) Earth model, which includes the effects of stratified structure, gravity, and Earth curvature. We note that whether the discontinuities must be included or not depends on the actual Earth model under consideration, e.g., model 1066A does not include discontinuities while 1066B does. In the case of a stratified model with discontinuities, the contribution of the discontinuities must be included but their contribution cannot be isolated. Based on the dislocation theory of the spherical Earth model, we calculate the effects of the discontinuities. The Love numbers and Green's functions caused by four independent point sources (strike-slip, dip-slip, horizontal tensile, and vertical tensile sources) are calculated using the 1066A and 1066B models. For the point sources, they are defined by UdS/R 2 = 1. Considering the discontinuity at 420 km, the Love numbers ( ) related to the sources at  419 km and 421 km are given in Figure 2. We need to calculate, to degree n, the maximum spherical harmonic degrees (n = 10.0·radius/source depth), to ensure sufficient precision. The 1066A model is continuous at 420 km, and the Love numbers of the sources at 419 km and 421 km are almost the same. However, for the discontinuous model 1066B, the Love numbers of the sources at 419 km and 421 km are quite different because the medium at 420 km is subject to changes due to the characteristics of a discontinuous interface.
h n , l n , k n l t n In Figure 2, are the Love numbers of spheroidal deformations, which represent the components of radial displacement, horizontal displacement, and gravitational potential, respectively.
is the horizontal displacement of toroidal deformation. The differences caused by the discontinuities (the blue lines) are large for the Love numbers of strike-slip and horizontal tensile sources, but for dip-slip and vertical tensile sources the differences are relat-l t n ively small. When the Love numbers for 1066A and 1066B are evaluated at radial distances above and below the discontinuity of 1066B, the differences obtained from 1066B are as expected much larger than those from model 1066A. The horizontal displacement of toroidal deformation , caused by a strike-slip source, shows a significant difference (almost 10%) in behavior between 1066A and 1066B. Due to the focal mechanism definition, toroidal deformations of tensile sources do not exist. In Figure 3, we show differences in synthetic Green's functions caused by four independent point sources at 419 km and 421 km, to evaluate the effects of discontinuities. is the Green's function for radial displacement, is the Green's function for horizontal displacement, is the Green's function for potential, and is the Green's function for gravity change at the surface. As expected, the difference between the Green's functions at 419 km and 421 km based on 1066B (the blue lines) is larger than that based . is the epicentral angle distance. is the radial displacement (row 1), is the horizontal displacement (row 2), is the potential (row 3), and is the gravity change at the surface (row 4).

Earth and Planetary Physics
doi: 10.26464/epp2021010 93 on 1066A (the red lines). For strike-slip and horizontal tensile sources, the difference based on 1066A is almost zero. For the dipslip source, the difference at epicentral angular distances of 0.02°< θ < 30° based on 1066A is not zero, but it is at least 10% smaller than the difference of 1066B. But for the radial displacement and gravity change caused by the vertical tensile source, the differences based on 1066A are three times as large as those based on 1066B at 0°< θ < 0.1°. In general and from all the components of the Green's functions, we find that the effects of the discontinuities are large.
Usually, a finite fault model can be divided into many subfaults. When the subfaults are defined small enough, they can be approximated as point sources Atkinson, 1997, 1998) and the individual contributions can be summed to get the total coseismic deformation of the finite fault. Any finite fault can thus be combined with the four independent point sources and its strike angle, dip angle, rake angle, etc. In view of the large effects of discontinuities on the Love numbers and Green's functions caused by point sources, we study below differences between application of the 1066A and 1066B models to the seismic deformations of the 2013 Okhotsk earthquake (M w 8.3), a deep-focus earthquake.

Effects of Discontinuities on Coseismic Deformations for the 2013 Okhotsk Earthquake (M w 8.3)
The main difference between the 1066A and 1066B models lies in their treatments of discontinuities at depths of 11.0 km, 21.0 km, 420 km and 671 km. To study the influence of the discontinuities in the 1066 model on seismic deformation, we take the 2013 Okhotsk earthquake (M w 8.3) as an example. The earthquake occurred on May 24, 2013, with its epicenter located at a depth of 598.1 km, below the Okhotsk Sea of Russia at (54.892°N, 153.221°E) (Figure 4a), between two discontinuous interfaces (420 km and 671 km). Although it caused little damage at the surface, the earthquake was felt across Asia as far as Moscow and across the Pacific along the western seaboard of the United States (USGS, 2013). Records of deep-focus earthquakes are rare in the world's modern seismograms, but their contribution to research on the effect of discontinuities is very prominent.
We use the slip model result inverted with a finite fault inverse algorithm (Ji C et al., 2002) by the United States Geological Survey (USGS, 2015) ( Figure 4b). Based on this slip model, we use the 1066A and 1066B models to compute the coseismic displacements and gravity changes caused by the Okhotsk earthquake. The calculated coseismic vertical displacements on the Earth surface using 1066A and 1066B models are shown in Figure 5(a, b). Figure 5c illustrates the difference between these two results. Figure 5 shows that the effect of discontinuities on the vertical displacements caused by Okhotsk earthquake appears mainly near the epicenter. The vertical displacements are negative east of the epicenter and positive west of the epicenter. However, the maximum and minimum values emerging from models 1066A and 1066B differ from each other; the values for the former are 9.49 mm and −16.4 mm, and the values for the latter are 9.72 mm and −15.5 mm, respectively. The differences around the epicenter are positive, reaching 2.1 mm. To give a relative numerical discrepancy ( ) caused by the discontinuities, we use Formula (1) to calculate the difference between the results of 1066A and 1066B. At the same time, we use a statistical method -root mean square (RMS) -for quantitative comparison, such as Formula (2). The discrepancy in vertical displacement due to the effects of discontinuities is 10.52%, with RMS error of 0.624 mm. The numerical results are summarized in Table 1.
Sea of Okhotsk where and are the computed coseismic surface deformations based on the 1066A and 1066B models, respectively, is the number of computed points, is the colatitude.
where and are the coseismic surface deformations of 1066A and 1066B, is the number of computed points.
The coseismic gravity changes for the 1066A and 1066B models are shown in Figure 5(d, e). Figure 5f displays the difference between these two models. The gravity decreases on the western side of the epicenter and increases on the opposite side. The maximum gravity changes are 3.01 μGal and 2.04 μGal for 1066A and 1066B, respectively. The minimum values are −1.57 μGal and −1.54 μGal, respectively. The differences around the epicenter are negative, reaching −0.17 μGal. Such small gravity change is expected to be detected by satellite mission GRACE (Xu CY et al., 2017). As shown in Table 1, the discrepancy in the gravity change is 6.19%, with RMS error 0.063 µGal. The geoid changes computed by models 1066A and 1066B (not shown here) exhibit similar differences. The differences around the epicenter are positive, reaching 0.072 mm. The discrepancy due to the effects of discontinuities is 9.07%, with RMS error 0.029 mm.
To further isolate the effects of how discontinuities on the far-field coseismic deformation are treated in Earth models, we present the far-field horizontal and vertical displacements as computed by the 1066A (blue arrows) and 1066B (red arrows) models and compare them with global positioning system (GPS) observations (yellow arrows) (Shestakov et al., 2014) in Figure 6. As for the 'far-field' definition, there is no unique standard. Usually, 'far-field' depends on fault size and seismic magnitude. We infer that it can be defined as distances from the seismic source that are at least equal to about twice the fault length. Due to the finite fault model inverted from GSN broadband data, not constrained by the GPS   data, the theoretical results are not in complete agreement with GPS data. Nevertheless, the theoretical displacement is generally consistent with the GPS data in both horizontal and vertical directions. The difference between the displacements of 1066A (blue arrows) and 1066B (red arrows) is quite small, with the two results almost overlapping.
The numerical RMS errors of far-field horizontal and vertical displacements between GPS and theories are presented in Table 2.
For the comparison between theoretical calculation (1066A) and GPS, the RMS error of horizontal displacements at N-S direction is 0.170 cm, while for E-W it is 0.252 cm. For vertical displacements, the RMS is 0.328 cm. Regarding the comparison between 1066B and GPS, the RMS of these three components is almost the same. The difference between 1066A and 1066B reflects the effects of the discontinuities, which are 0.02 cm, 0.037 cm and 0.032 cm, respectively. The difference between the 1066A and 1066B models is small, compared with the GPS displacements. Consequently, we   can neglect the effect of discontinuities on far-field coseismic deformation but not on near-field deformation.

Effects of Discontinuities on Coseismic Deformations for the 2011 Tohoku-Oki Earthquake (M w 9.0)
Since most earthquakes occur in the shallow crust, i.e. the top 100 km, accurate calculation of Green's functions caused by sources within this depth range is necessary. We use models 1066A and 1066B to study the coseismic deformations of the 2011 Tohoku-Oki earthquake (M w 9.0), a shallow event. The geodetic observations, GPS, and Gravity Recovery and Climate Experiment (GRACE), can detect gravity changes due to mass rearrangements caused by the earthquake. Comparing such data with the computational results from dislocation theory, we can quantitatively estimate the effects of the discontinuities over this class of data.
The 2011 Tohoku-Oki earthquake (M w 9.0) occurred at a depth of 24.4 km near the east coast of Honshu in Japan on March 11, 2011, with its hypocenter located at (38.322°N, 142.369°E). Japan is one of the world's most earthquake-prone countries, situated atop four tectonic plates: the Eurasian, the Pacific, the Okhotsk, and the Philippine Plates. This 2011 event is the most devastating earthquake in the modern seismic record of Japan. Following the earthquake, Japan experienced a very large tsunami, destructive fires, and destruction of nuclear power facilities. We take the fault slip model from Advanced Rapid Imaging and Analysis (ARIA) (Wei SJ et al., 2011), inverted with the GSN broadband data and constrained by GPS data, to calculate the coseismic deformations. In order to estimate the effects of the discontinuities, we make use of the slip model from the USGS website, inverted from GSN broadband data only.
After the 2011 Tohoku-Oki (M w 9.0) earthquake, the data from the GPS Earth Observation Network (GEONET) covering all of Japan showed that the island moved eastward, with a maximum displacement of 5.3 m. GPS results from the Crustal Movement Observation Network of China (CMONOC) also detected far-field deformation (Wang M et al., 2011) in eastern China. Since the effects of discontinuities are negligible on far-field deformation, we focus on near-field deformation.
We use the 1066A and 10666B models to compute the near-field coseismic vertical and horizontal displacements. The theoretical calculation and the GPS data are presented in Figures 7 and 8. The displacements calculated with the ARIA slip model are closer to the GPS data than those from the USGS model. The consistency is largely due to the fact that the GPS data are used to invert the ARIA slip model. The horizontal displacements based on model 1066B are in closer agreement with the GPS data, opposite to the vertical displacements from the statistical analysis (Table 3). We calculate the RMS error to make the statistical analysis. The RMS caused by the discontinuities on horizontal displacements are 0.030 m (N-S) and 0.093 m (E-W); on vertical displacements the RMS is 0.025 m (up-down) in the study area computed with ARIA slip model, third row of Table 3. For the results from the USGS model, the effects on horizontal displacements are 0.011 m (N-S) and 0.026 m (E-W) and 0.013 m on vertical displacements (updown), respectively, sixth row of Table 3. Note that, in the inver-sion of the slip, if the Earth model is the same as that used in the calculation to evaluate the deformation and contains the discontinuities, the inverted slip model is expected to be more precise when fitting the GPS data.
We also use the ARIA and USGS slip models to compute the vertical displacements and gravity changes (Figure 9), the latter being more closely linked to the vertical than to the horizontal. Figures  9(a, b) and 9(g, h) correspond to the vertical displacements on the Earth surface calculated by the 1066A and 1066B models, respectively. Figures 9(c, i) provide differences between the two models. Figures 9(d, e) and Figure 9(j, k) are the gravity changes for the 1066A and 1066B models. Figures 9(f, l) are the differences. The tendencies of variations in vertical displacements and gravity changes for these two models are consistent but differ in terms of the maximum and minimum values. Both models show positive changes west of the epicenter and negative changes east of the epicenter for gravity changes, opposite to the vertical displacements. The statistics of Figure 9 are shown in Table 4. The discrepancy in the vertical displacement, which reflects the difference caused by the discontinuities, is remarkably large, reaching 12.1% for the ARIA model and 11.6% for the USGS model. The discrepancy in the gravity change is 11.8% for the ARIA model and 10.7% for the USGS model. For the geoid changes (not shown here), it is 10.4% for the ARIA model and 7.5% for the USGS model.
To investigate the effects of discontinuities on the gravity change in the spatial domain we calculate the coseismic gravity changes at spatially fixed points. We divide the study area into grids of 1°×1°, comparable to the resolution of the GRACE measurements.
In this study we use monthly CSR GRACE Release 05 (RL05) Level 2 data products (Save et al., 2016) from January 2003 to May 2011 to calculate the gravity changes caused by the 2011 Tohoku-Oki earthquake (M w 9.0). We take 95 months of data, from January 2003 to December 2010, to determine the gravity field.
Then, with the 60 degrees Stokes coefficients from GRACE, we calculate the monthly gravity changes and remove the mean-field gravity (from January 2003 to December 2010) at spatially fixed points for the area from 130°E−155°E longitude and 25°N−50°N latitude with grids of 1°×1°. We take the mean values for March to June annually (2003−2010) to represent the pre-seismic deforma-   al 'stripes', we apply a filtering scheme with decorrelation filtering P3M6 (a third-degree polynomial fitting for order 6 and above of Stokes coefficients) (Swenson and Wahr, 2006) and 300 km Gauss smoothing (Wahr et al., 2004;Han et al., 2005). Figure 10a shows the gravity changes caused by the earthquake. The variation characteristics are similar to the results of Matsuo and Heki (2011), but the minima are different because these authors used the GRACE RL4 data and a different data processing approach. The largest To compare the theoretical results obtained from Earth models with the GRACE signals, we compute the gravity changes at spatially fixed points and truncate them to the same degree as GRACE, up to degree and order 60, and apply decorrelation filtering P3M6 and 300 km Fan filtering (Zhang ZZ et al., 2009). Here we compute the gravity changes based on the ARIA slip model, whose results are closer to the observation data. Because most crustal uplift occurred in the ocean, we also make the seawater correction (De Linage et al., 2009;Broerse et al., 2011;Sun WK and Zhou X, 2012), using the Spherical Bouguer correction (Yang JY et al., 2015) to compute the effects of seawater mass. The results are presented in Figure 10. Figure 10(b, c) shows the calculated coseismic gravity changes (up to degree and order 60) obtained from models 1066A and 1066B. Figure 10d displays difference between results of the two models. Figure 10(e, f) presents differences between GRACE data and modeled data (1066A, 1066B), respectively. The main difference between the two theoretical models' results occurs in the Pacific Ocean. For the statistical results, the RMS error between observed (GRACE) and modeled 1066A is 1.378 µGal and for model 1066B is 1.377 µGal. The RMS error caused by the discontinuities is 0.0016 µGal and the discrepancy is 11%. The effects of discontinuities cannot be seen at the spatial resolution of GRACE data. But the theoretical calculation from 1066B is closer to GRACE than that from 1066A.

Discussion
The philosophy underlying this study is that all investigators should be aware of the potential effects of how discontinuities are treated in various Earth models, because different approaches can lead to significant differences in computed coseismic displacements and gravity changes.
Accordingly, this study has investigated how differences between how SNREI Earth models treat discontinuous Earth structures tend to affect computed coseismic displacements and gravity changes.
Our work is a theoretical study to call attention to the importance of the effects of Earth's discontinuities in Love number and Green' s functions evaluation, specifically within the framework of the 1066 model, comparing 1066A with 1066B.
The Love numbers and Green's functions from four independent point sources are computed for the 1066A and 1066B Earth models. The effects of the discontinuities on Love numbers of strikeslip and horizontal tensile sources are found to be large, compared to cases of dip-slip and vertical tensile sources. The difference between models 1066A and 1066B in computed horizontal displacement in toroidal deformation, caused by a strike-slip source, is almost 10%. For most components of Green's functions, the difference between the sources at 419 km and 421 km based on 1066B is larger than 1066A. We thus conclude that the effect of discontinuities is worth taking into consideration when applying 1066 Earth models.
Earth modeling is a dynamic field of study. The latest PREM (Dziewonski and Anderson, 1981) and AK135 (Kennett et al., 1995) models have become widely used, approximating well the real Earth's radial stratification. For example, PREM consists of 96 layers, including 12 discontinuous interfaces. But the increase in the number of discontinuous interfaces has led researchers to reduce calculation time by using numerical interpolation methods in the dislocation theory. The discontinuities may thus be ignored or smoothed through interpolation, but our results clearly indicate that discontinuities inside Earth models should be considered self-consistently.
It is worth noting that this study has considered the layered spherically symmetric Earth model. A three-dimensional Earth model would be more realistic, but 3D Earth models are more complex.
Scientists usually use the finite element method to simulate and calculate seismic deformation of the 3D Earth model. On the other hand, the seismic deformation theory of the three-dimensional Earth model, although very complex, can also be used (for example, Fu GY and Sun WK, 2007). In this case, the lateral inhomogeneity effect and the interfault effect will be coupled together, but we can still calculate and discuss them. In fact, to do so we need only to add three-dimensional structural parameters to the present layered Earth model (1066 or PREM).

Conclusions
The coseismic deformations for the 2013 Okhotsk earthquake (M w 8.3) indicate that the 1066A and 1066B models must be investigated regarding their different behavior in terms of the effects of Earth's discontinuities. The models' discrepancies (between calculated effects of discontinuities) in coseismic vertical displacement, geoid, and gravity changes are 10.52%, 9.07%, and 6.19%, with    These results indicate that the effects of discontinuities on the near-field coseismic deformations are large and nonnegligible. For the far-field deformation, differences between models 1066A and 1066B in their treatment of discontinuities lead to differences in calculated horizontal displacements (N-S, E-W), and vertical displacements (up-down) of 0.02 cm, 0.037 cm, and 0.032 cm, respectively, which are small and negligible.
For a shallow event, we choose the 2011 Tohoku-Oki earthquake (M w 9.0), which occurred within 100 km depth. The near-field displacements evaluated via the 1066B model are in better agreement with the GPS data than those from 1066A. In addition, the displacements calculated with ARIA slip model are closer to the GPS than those from USGS. For the results from the ARIA slip model, the effects in our study area of the discontinuities are: on horizontal displacements, 0.030 m (N-S) and 0.093 m (E-W); on vertical displacements, 0.025 m (up-down). The discrepancy (the effects of discontinuities) is approximately 12.1%, 10.4% and 11.8% for vertical displacements, geoid, and gravity changes, respectively. The GRACE-detected gravity changes for this earthquake are −6.2 μGal in the subduction zone and 2.3 μGal in the ocean. The RMS error between observed (GRACE) and modeled 1066A is 1.378 µGal; for modeled 1066B, it is 1.377 µGal. The difference due to the effect of discontinuities is 0.0016 μGal and the discrepancy is 11%. The modeled result from 1066B is closer to agreement with empirical GRACE data.   In summary, the effects of discontinuities on coseismic deformations are quite large for point sources and finite fault models. Our results suggest that, for most earthquakes, the treatment of discontinuities inside an Earth model should be carefully addressed when computing coseismic gravity and displacements within spherical dislocation theories, especially for the near-field; Green's functions on both sides of the discontinuities should be computed separately regardless of the numerical scheme used.