Co-seismic slip distribution of the 2011 Tohoku (M_W 9.0) earthquake inverted from GPS and space-borne gravimetric data

1.   Introduction

Earthquakes radiate seismic waves and cause surface displacements and gravity anomalies associated with the redistribution of Earth mass. Thanks to the rapid development of space geodetic techniques, gravity anomalies caused by megathrust earthquakes have been detected by the GRACE (Gravity Recovery and Climate Experiment) satellite mission. Three of the most important seismic events documented by GRACE include the 2004 Sumatra M_W9.3 earthquake, the 2010 Chile M_W 8.8 earthquake, and the 2011 Tohoku M_W 9.0 earthquake (e.g., Han SC et al., 2006; Chen JL et al., 2007; de Linage et al., 2009; Heki and Matsuo, 2010; Zhou X et al., 2011; Matsuo and Heki, 2011; Zhou X et al., 2012). The GRACE satellite mission provides time-series data describing the long-wavelength gravity field over a time period longer than 10 years. Theoretically, it can detect static co-seismic gravity changes induced by events greater than 8 in magnitude (Gross and Chao BF, 2001; Sun WK and Okubo, 2004). However, due to noise sources peculiar to GRACE data, no earthquake with a magnitude less than 8.8 has yet been detected.

Based on the forward elastic dislocation model, co-seismic changes, including changes in geoid shape, deflection, and gravity gradients observed by GRACE, have been modeled using an a priori focal mechanism inferred from seismological or geodetic data (e.g., de Linage et al., 2009; Sun WK and Zhou X, 2012; Zhou X et al., 2012; Wang L et al., 2012a). Han SC et al. (2011) constrained the dip angle, seismic moment, and depth of the 2011 Tohoku earthquake using inter-satellite range rate measurements made before and after the event. Wang L et al. (2012b) implemented the orthogonal Slepian functions to extract from GRACE time-series data the co-seismic gravity changes associated with the 2010 Chile M_W 8.8 event. They then inverted the fault plane size and the uniform slip using the simulated annealing method. Cambiotti and Sabadini (2012) and Wang L et al. (2012c) applied similar methods to the 2011 Tohoku earthquake. Later, Cambiotti and Sabadini (2013) obtained the first Gravitational Centroid Moment Tensor (GCMT) solution of the 2011 Tohoku earthquake. This solution, based solely on GRACE data, consists of the centroid and the seismic moment tensor, as well as the associated error bounds. Dai CL et al. (2014) used the change in the north component of gravity and corresponding gravity gradients to obtain both the rake angle and the centroid location. In this regard, GRACE improves our understanding of the physics of the earthquake source. Most recently, Fuchs et al. (2016) investigated the fault slip distribution for the 2011 Tohoku earthquake using data obtained using GPS, GRACE, and GOCE (the Gravity Field and Steady-State Circulation Explorer).

Various geophysical measurements, including seismic waveforms, ground motions, displacements, and tsunami events, have been used to characterize finite fault models of megathrust earthquakes with heterogeneous slip distributions. Such models, which describe the slip distribution over the fault in detail, form the basis of investigations into stress changes, which have utility in the realms of seismic hazard mitigation and earthquake cycle modeling. In the case of the 2011 Tohoku earthquake, the largest inter-plate event on the Pacific-Okhotsk convergent boundary in several centuries, many finite fault models have been obtained from seismic waveforms (e.g., Lay et al., 2011; Hayes, 2011; Ide et al., 2011; Lee et al., 2011), geodetic displacements (e.g., Ozawa et al., 2011; Pollitz et al., 2011; Loveless and Meade, 2011), and combined seismic waveforms, geodetic data, and tsunami observations (e.g., Simons et al., 2011; Yagi and Fukahata, 2011; Yokota et al., 2011; Koketsu et al., 2011; see more source inversion models in the paper of Tajima et al., 2013). Although large differences in the slip distributions exist among these models due to the differing methods for dealing with the different data sources, data describing unique seafloor deformations measured by GPS-Acoustic (GPS-A) (Sato et al., 2011; Kido et al., 2011) were used to improve the slip resolution in the shallow part of the fault. The results indicate that up-dip slip exceeding 40 m occurred near the trench (Zhou X et al., 2014).

In this paper, we extend the inversion of GRACE gravity data to the case of heterogeneous slip during the 2011 Tohoku earthquake. We discuss the resolving power of GRACE on the fault plane and compare, within the same inversion scheme, results obtained using GRACE data to results obtained using other types of geodetic data, such as GPS displacements. Furthermore, we discuss the effects of fault geometry on the earthquake gravity pattern and the resulting slip inversion. Finally, we jointly invert from GPS and GRACE data to obtain a slip model.

2.   GPS and GRACE Data

2.1   GPS Data

To better understand the resolving power of GRACE gravity data, we herein compare it to continuous GPS (cGPS) co-seismic displacements and jointly invert for the slip distribution over a finite fault using both displacement data and gravity changes. We consider the GPS co-seismic displacement (version 0.3) determined by ARIA at the Jet Propulsion Laboratory based on GEONET data provided by the Geospatial Information Authority (GSI). The co-seismic jumps of 1232 cGPS stations employed in Japan were determined using offset data recorded 6 minutes before the main shock and 9 minutes after. These data, which include no aftershock-related displacements, have been used by many researchers (e.g., Pollitz et al., 2011; Zhou X et al., 2014). The onshore GPS data describe a maximum eastward displacement of 5 m. Moreover, data from 7 GPS-Acoustic (GPS-A) stations, describing a maximum horizontal displacement of 31 m, improve the spatial resolution on the shallow part of the fault. These data are used for joint inversion (Sato et al., 2011; Kido et al., 2011). It is important to note that the GPS-A data have large uncertainties, of 0.5 m and 1.0 m for the horizontal and vertical components, respectively.

2.2   Co-seismic Gravity Changes by GRACE Measurements

To estimate the co-seismic gravity change associated with the 2011 Tohoku earthquakes, we used GRACE Level2 RL05 time-series data recorded between January 2003 and December 2013 and released by the Center for Space Research (CSR). They are given in terms of Stokes coefficients up to harmonic degree 60 and with a temporal resolution of one month. For our analysis, data for Stokes coefficient C₂₀, related to the Earth’s oblateness, are replaced by equivalent data measured using Satellite Laser Ranging (SLR) due to the poor accuracy of GRACE mission data for this coefficient (Cheng MK and Tapley, 2004). The gravity signals of ocean tides, pole tides, and atmospheric effects have all been removed using GRACE RL05 products. The hydrological signal was not removed. Although the latter may be significant, tests performed using the GLDAS hydrological model (Global Land Data Assimilation System (Rodell et al., 2004)) showed negligible effects of the hydrologic signal on the estimate of the co-seismic gravity change. Since the highest hydrological gravity effect in the study is only 0.5 μGal (over the hanging wall region, where the largest gravity signal is attained, as shown in Figure S1), we decided not to remove the signal. The largest hydrological effect is attained on the northwestern boundary of the study area, away from the region affected by the earthquake.

When using GRACE temporal gravity field data, it is common to reduce the significant stripes resulting from high-frequency noise using a post-processing scheme based on the combination of empirical polynomial fitting decorrelation and a Gaussian filter (Swenson and Wahr, 2006). Despite its popularity, this scheme has a remarkable disadvantage: the polynomial is selected artificially, and it might dampen the seismic gravity signals and change the spatial pattern orientation (Dai CL et al., 2014). To overcome the above shortcoming and to enhance the signal-to-noise ratio, we use a non-isotropic DDK3 decorrelation filter (Kusche, 2007; Kusche et al., 2009).

We calculate the time series of gravity changes using spherical harmonic summation in the spatial domain over a 0.5° grid for the study region, which is a rectangular area reaching from 130°E to 155°E and from 25°N to 50°N. At each point of the grid, we fit the gravity time-series data to a linear model (Appendix A), which includes long-term trends, annual, semi-annual and 161 day S2 ocean tidal aliasing periods, a Heaviside function that describes the co-seismic gravity change, and an exponential function that quantifies the post-seismic gravity decay (See Figure 1a, 1b). The resulting co-seismic gravity change map and its error distribution are shown in Figures 1(c, d). The co-seismic gravity change varies from -10.6 to 3.6 μGal over the selected grid points and is characterized by an asymmetric dipolar pattern toward the negative pole (Figure 1c).

Besides the long-term, periodic, and co-seismic signals, GRACE time-series data also include a post-seismic signal. This signal can be seen as a drop or a jump in the gravity change in the first three months after the quake, which rebounds as an exponential curve in the following months (Figures 1(a, b)). A similar pattern in time was found for the 2004 Sumatra and the 2010 Chile earthquakes by Tanaka and Heki (2014), who used two exponential functions to fit the data. The physical mechanism of this signal is not clear and to model it would be beyond the scope of this study. Despite this, we must remove the post-seismic signal from the time-series because we have found that it can contaminate the co-seismic jump and affect the resulting fault slip inversion. This finding differs from those of Dai CL et al. (2014), who claimed that the post-seismic signal has only a minor impact on estimates of co-seismic change. In time-series analysis, we use an exponential function with a uniform relaxation time of 10 months over the map to describe the feature of post-seismic decay by a considerable number of tests. Taking the great disturbance in the month of event occurrence into account, we ignore the month of March 2011.

After spatial smoothing, we initially consider an a priori uncertainty of 1 μGal for each selected grid. Using the formula (A1), the post-fit uncertainties of the co-seismic gravity change-fitting model are estimated in terms of the rescaling factor ${\chi ^2}$ . The posterior uncertainties that do not exceed 2.5 μGal are then used to define the weight matrix for the slip inversion (Figure 1d). The largest errors are in areas further than 45° north of the equator. In this study, the dipole-patterned gravity changes in the epicenter area, where the signal-to-noise ratio (SNR) is highest, will be used to constrain the fault slip pattern.

Figure  1.  Panels (a) and (b) show the time-series analysis of GRACE data for the 2011 Tohoku earthquake at points A (139°E, 39°N) and B (143°E, 34°N), which are denoted by white triangles in panel (c). The original monthly gravity change is denoted by the black line and associated black points, while the linear fitting model is denoted by the blue lines (which include the long-term trend, periodical signals, the co-seismic jump, and the post-seismic exponential decay) and the red lines (which include only the long-term trend, the co-seismic jump, and the post-seismic exponential decay). Panels (c) and (d) show the co-seismic gravity changes estimated from GRACE data analysis and their post-fit uncertainties in the region under study. The black star indicates the epicenter of the earthquake, which was obtained from the Japan Meteorological Agency (JMA).

DownLoad: Full-Size Img PowerPoint

3.   Method

3.1   Finite Fault

To obtain a static slip model of the 2011 Tohoku earthquake, we must first discretize a finite fault. Previous work has shown that the subduction angle of the Pacific plate at the Japan Trench becomes steeper with depth (Tsuru et al., 2002; Takahashi et al., 2004), and that the trench-dependent strike angle varies between 190° and 210° (Hayes et al., 2012). However, results obtained using waveform and GPS data show that slip inversions do not differ markedly between a simple planar fault and a complex non-planar fault (Lee et al., 2011; Zhou X et al., 2014).

Thus we first consider a simple rectangular fault plane (I) with an along-strike length of 750 km and an along-dip length of 252 km. This fault is divided into 30 × 14 patches of 25 km × 18 km each. The fault plane geometry, as shown in Figure 2a, is obtained by averaging the Slab1.0 subduction model (Hayes et al., 2012); the resulting fault is characterized by dip and strike angles of 11° and 201°, respectively. To analyze the sensitivity to fault geometry of slip inversion results based on GRACE data, we will consider a second fault (II), similar to Fault I but with a variable dip and realistic non-planar geometry, shown in Figures 2b and 2c, respectively. The dip of Fault II varies from 3° to 22° and tends to increase with increasing depth. Faults I and II are similar to those described by Zhou X et al. (2014) but with a larger fault size along the trench.

In fact, the strike along the trench is variable. To check the sensitivity of gravity data to the strike, we construct a non-planar fault with a real strike. On the basis of Slab 1.0, the strike of Fault III varies with trench orientation from 190°–210°, and the dip varies with depth between 3° and 25° (Figure 2c). The dip angle is like that of Fault II; however, the width of the fault is extended to 320 km. We divide Fault III into 30 × 16 patches of 25 km × 20 km. Relative to the strikes of Faults I and II, the variable strike of Fault III is more coincident with the trench; this geometry is more realistic for the subduction zone.

To overcome the computational difficulty of implementing Green’s function evaluation within a spherically layered model, we set the top edge of all the three faults to 2 km below the seafloor. Note that the dimension of each patch should not be confused with the resolving powers of GRACE or GPS since these are defined by the interference of all the patches collectively, as we show in the section 4.1.

Figure  2.  The three finite faults used to invert the slip distribution. (a) Fault I is a simple fault with a uniform strike of 201° and a dip of 11°. The fault boundary and subfaults are denoted by black rectangles and red points, respectively. Fault II (b) is similar to Fault I but with a variable dip ranging from 3° to 22°. Fault III (c) is characterized by almost the same area as Fault I. However, Fault III is based on more realistic geometric parameters, including variable strike and dip, obtained using the Slab 1.0 subduction model.

DownLoad: Full-Size Img PowerPoint

3.2   Inversion Method

For both the along-strike and dip-slip components of each patch, we derive the Green’s function, which describes co-seismic displacements and gravity changes on the same grid that was used for the GRACE data analysis (after DDK filtering). These Green’s functions are obtained within the framework of a spherically symmetric, self-gravitating compressible Earth model (Sun WK and Okubo, 1993; Sun WK et al., 2009) based on the Preliminary Reference Earth Model (Dziewonski and Anderson, 1981). They account for the direct contributions of ocean water redistribution to changes in gravity by means of an approximate infinite plane (Bouguer layer) model (see Appendix B for details). In this respect, we neglect the loading of the ocean water redistribution because this second-order effect is small (Broerse et al., 2011).

Because they take self-gravitation and the radial stratification of the material parameters of the Earth into account, our Green’s functions are more realistic than traditional functions obtained within the framework of homogeneous half-space Earth models (Okubo, 1992; Han SC et al., 2006; Wang L et al., 2012b, c).

For consistency with the GRACE scheme, we use the spherically layered dislocation based on PREM to compute the displacement Green’s function for the GPS data.

The inversion method for obtaining the slip distribution minimizes the cost function

where d is the observation data, G is the Green’s function, s is the slip distribution, ${{\varSigma }}$ is the covariance matrix, L is the regularization matrix based on the second order Laplacian operator (Harris and Segall, 1987), and $\alpha$ is the regularization parameter, which we determine by finding the knee of the L-curve (Hansen and O'Leary, 1993). A non-negative constraint is used to avoid unreasonable rake angles (Lawson and Hanson, 1974). We adopt the inversion scheme based on rigid boundary conditions, meaning that large slips do not reach the trench but that slip does occur on the shallow part of the fault (Zhou X et al., 2014).

4.   Results

4.1   Checkerboard Tests

To better understand the resolving powers of GPS and GRACE, respectively, before the slip inversion using real data, we check the resolutions of the GPS displacements and the long-wavelength (DDK3-filtered) GRACE gravity data on the fault plane using the checkerboard test. The prior checkerboard model, with an alternating up-dip slip of 20 m and 0 m over Fault I, has a slightly larger seismic moment than that of the 2011 Tohoku earthquake. The results shown in Figure 3 indicate that the resolving power of co-seismic gravity data is about 300 km (corresponding to about half the wavelength of the GRACE data) over the entire fault depth, including the shallow part of the fault, where large amounts of slip did indeed occur. However, at the deepest parts of the fault inland, cGPS performed better (Zhou X et al., 2014). Once we combine the synthetic cGPS and the synthetic GRACE data to invert the fault slip, the inverted slip resembles that obtained using synthetic cGPS data. Thus, adding cGPS data to GRACE data can be expected to improve the spatial resolution when we resolve the slip distribution of a megathrust event using space-borne gravimetric observations. In Figure 3e, we combine GRACE, cGPS and GPS-A data, all of which are synthesized by the input slip model, to test joint inversion performance. Because of the improved resolving power of GPS-A on the shallow part of the fault, the jointly inverted slip more closely matches the slip pattern near the trench than the slip calculated using data without GPS-A measurements (this can be seen by comparing the results for 5×2 subfaults. This finding matches previously documented findings (e.g., Ozawa et al., 2012; Zhou X et al., 2014).

Note that the weight ratio among these data impacts the slip inversion. Here, we adopt the same weight of 1 : 10 : 10 for cGPS, GPS-A and GRACE, as is used in the real data inversion that follows. These findings thus suggest that in order to improve the reliability of slip inversion throughout the whole fault depth by using GRACE data and enhance the spatial resolution, we must use cGPS or/and GPS-A data on the solid earth surface.

Figure  3.  Checkerboard test. Left panels (a) show the prior slip models with alternating up-dip slips of 20 m and 0 m (black and white blocks) over 2×2, 3×2 and 5×2 subfaults (top to bottom). Panels (b) and (c) show the inverted slip model from the synthetic displacements and co-seismic gravity changes obtained using the prior model over 2×2, 3×2 and 5×2 subfaults. Panels (d) and (e) show the joint inversion results obtained using synthetic GRACE and cGPS datasets and GRACE, cGPS, and GPS-A datasets, respectively. The Green’s function for gravity changes and the added Gaussian noise are smoothed by means of the DDK3 filter.

DownLoad: Full-Size Img PowerPoint

4.2   Slip Inversion from GPS

In order to analyze the concordance between the slip patterns inverted from GRACE and those inverted from GPS, we invert for the finite slip over the three fault models (slip-cGPS-I, slip-cGPS-II and slip-cGPS-III) from cGPS data with a greater resolution following the scheme described in Zhou X et al. (2014). As shown in Figure 4, although Faults I and II are extended in along-strike length compared to those described in Zhou X et al. (2014), the slip patterns are fully consistent with those described in that study. As pointed out by Zhou X et al. (2014), cGPS is only slightly sensitive to the dip of the subducting plate. Compared with slip model results obtained for Faults I and II, results obtained for Fault III (Figure 4) exhibit an important difference: the slip distribution is more localized (close to the epicenter). This can be seen in the overlap of the highest slip values with the black star (bottom left panel of Figure 4 compared). The seismic moments for the three cGPS slip models are 3.55×10²² Nm, 3.54×10²² Nm and 3.00×10²² Nm, respectively (Table 1). Each slip model containing values exceeding 30 m has a significant fan-like rake pattern. The slip models agree with those developed in other cGPS studies (Ozawa et al., 2011; Pollitz et al., 2011; Loveless and Meade, 2011; Zhou X et al., 2014).

Using the slip models inverted from cGPS, we can synthesize the spatial gravity changes with the Green’s function described in section 3.2. These results can be compared with GRACE measurements. The second column of Figure 4 shows that similar dipolar patterns of gravity changes in Figure 1 are predicted by slip-cGPS-I, slip-cGPS-II, and slip-cGPS-III. The patterns of the three modeled gravity changes resemble GRACE observations with correlation coefficients of 0.89, 0.90, and 0.91, respectively. Despite the slight improvement in the correlation coefficients, these cGPS slip models can explain only about 50% of the gravity amplitudes (third column of Figure 4). The large discrepancies between modeled and observed data can be ascribed to differences in temporal resolution: co-seismic jumps in GPS data are measured minutes before and after the main shock, while GRACE data are obtained in monthly time-steps.

Figure  4.  Co-seismic slip distributions on Faults I, II, and III (shown in left panels from top to bottom) are inverted from cGPS displacements. Each contour represents 5 m slip. Co-seismic gravity changes modeled by the slip models in each left panel are shown in the corresponding middle panels. The corresponding right panels show differences between model results and GRACE observations (Figure 1c).

DownLoad: Full-Size Img PowerPoint

4.3   Slip Inversion from GRACE

4.3.1   Slip model on Fault I

By the same method we used for the checkerboard test, we can obtain the slip distribution, shown in Figure 5a, for the planar Fault I (slip-grac-I) from inversion of the co-seismic gravity change obtained from GRACE. The total seismic moment is 4.9×10²² Nm (equivalent to M_w=9.06), the mean rake angle is 80.9°, and the maximum slip amplitude is 22 m. The pattern is smooth, and slip is concentrated in the central part of the fault. The centroid (represented by the empty triangle) is displaced slightly to the southwest of the JMA epicenter. The high-slip area in Figure 5a is also displaced southward, but it extends more widely along the strike than the cGPS slip pattern (Figure 4, top left panel). The area where the rupture exceeded 10 m has a length of about 500 km along the strike and a width of about 180 km along the dip. This is in agreement with the checkerboard test results, showing a uniform resolving power with depth; the slip is characterized by rake angles of about 90°. Slip amplitudes larger than 5 m are present over the entire fault.

The rake pattern obtained using GRACE data is closer to pure up-dip than the pattern obtained using inland GPS data, which is consistent with GPS-A observations (Zhou X et al., 2014). The modeled co-seismic gravity change, shown in Figure 5b, agrees in terms of orientation of the largest notch of positive gravity anomalies with the GRACE gravity change shown in Figure 1a. Gravity varies from -10.6 to 3.5 μGal and accounts for 70% of the co-seismic gravity change measured by GRACE; see Figure 1a. Small discrepancies of at most 6% of the largest peak-to-peak gravity anomaly of Figure 1 can be seen between modeled and observed gravity changes in the epicentral area. These discrepancies can be seen in the residuals of Figure 5c. Gravity changes in the surroundings of the dipolar pattern cannot be completely explained by slip-grac-I. The most significant discrepancy appears in East Asia, at the margin of the data sampling area. Evidently, this discrepancy is not related to the earthquake.

By inverting uniformly sampled GRACE gravity data, we obtain a physically sound slip distribution I on the Fault I plane for the 2011 Tohoku earthquake, portrayed in Figure 5a. The maximum slip is 22 m, and the mean rake angle is 81°; these values are consistent with a thrust earthquake, and the rake angles approach those inferred by other GRACE studies (shown in Table 1). The total seismic moment of slip-grac-I is smaller than the seismic moment estimated by other studies relying on GRACE data, and the centroid determined by slip-grac-I is located at a greater depth (20 km, Table 1). This discrepancy can be ascribed both to the trade-off between seismic moment and depth (Han SC et al., 2011) and to the fact that the results of the other studies consider a point-like seismic source or a fault with uniform slip. Our model, which considers a heterogeneous slip distribution, allows us to account for a more realistic seismic moment distribution with depth. In this respect, it is noteworthy that our estimate of the seismic moment of slip-grac-I is in closer agreement with estimates obtained from the inversion of other datasets, such as seismic waveforms (Hayes, 2011) and GPS (Ozawa et al., 2011; Zhou X et al., 2014), than estimates made in other studies based on GRACE data (Table 1). The fact that our estimate is larger than all the estimates obtained using non-GRACE datasets may be ascribed to the monthly temporal resolution of GRACE, which includes changes that may have occurred after the event (Han SC et al., 2014).

Our slip-grac-I data (Figure 5a) are characterized by a smooth and extended distribution with respect to data obtained by inverting inland GPS measurements (Figure 4, top left panel). Compared to results obtained by Dai CL et al. (2014), which inverted the north component of gravity changes and gravity gradient changes, our results for the centroid more closely match the actual epicenter as measured by the JMA (142.8°, 38.05°).

Figure  5.  Slip-grac-I (a) inverted from co-seismic gravity changes estimated from GRACE data (Figure 1c); each contour line represents 5 m slip. Co-seismic gravity change (b) modeled by the inverted slip distribution and associated residuals (c) with respect to co-seismic gravity changes as estimated from GRACE data. The triangle represents the centroid of the fault slip distribution, which is shifted to the southwest of the JMA epicenter (star).

DownLoad: Full-Size Img PowerPoint

The slip-grac-I shown in Figure 5a was inverted using the best smoothing parameter ( $\alpha$ = 0.501), which is located at the knee of the trade-off curve (Figure 6a). As shown in Figures 6b and 6c, rougher and smoother slip models are inverted for smoothing parameters $\alpha$ = 0.391 and $\alpha$ = 0.631, respectively. As expected on the basis of Figure 3, both portray the same power to resolve depth. We use the L-curve to select an optimal smoothing factor to obtain the slip distribution in the following inversion.

Figure  6.  (a) The regularization parameters used during slip model inversion shown in the trade-off curve. (b) Rougher and (c) smoother slip models inverted from GRACE-observed gravity changes using regularization parameters of 0.391 and 0.631, respectively. Each contour represents 5 m of slip.

DownLoad: Full-Size Img PowerPoint

4.3.2   Slip models on Fault II

Using the same scheme as slip model I, we invert a co-seismic slip distribution on Fault II (slip-grac-II) using GRACE gravity change data with a uniform spatial sampling scheme. The result is shown in Figure 7a. The total seismic moment of slip-grac-II is 5.1×10²² Nm (equivalent to Mw=9.07). This value is slightly larger than that obtained using slip-grac-I but smaller than those obtained in other studies relying on GRACE. The maximum slip is about 20 m, and the average rake is 80°, approximating Han SC et al. (2011) and Dai CL et al. (2014) (Table 1). The slip distribution on Fault II, which has a non-uniform dip, resembles the distribution found for slip-grac-I, shown in Figure 5a, although the distribution does trend slightly deeper due to the increased slab dip, and its centroid has moved slightly southward. The fan-like rake pattern is consistent with the patterns obtained using slip-grac-I and from GPS data. As was the case for results obtained using slip-grac-I, the residuals between modeled and observed gravity changes shown in Figure 7c are very small in the dipolar pattern region. Slip-grac-II can explain 70% of the GRACE measurements, the same as can slip-grac-I.

Comparison of slip-grac-I and slip-grac-II indicates that the fault dip angle affects only the seismic moment; it does not impact the slip pattern inverted from the GRACE gravity change data. In other words, results based on GRACE gravity observations are only slightly sensitive to the subduction angle of the Pacific slab. In the next section, by means of the slip distribution for a more realistic fault with realistic strike changes along the Japan Trench, we will verify whether the fault strike affects the slip inversion.

Figure  7.  Slip distribution on Fault II (Figure 2b) (slip-grac-II (a)) inverted from co-seismic gravity changes estimated from GRACE data (Figure 1c). Co-seismic gravity changes predicted by inverted slip model II and the associated residuals are shown in (b) and (c).

DownLoad: Full-Size Img PowerPoint

4.3.3   Slip models on Fault III

By the same approach as before, we obtain the co-seismic slip model on Fault III (shown in Figure 8 as slip-grac-III), which has a more realistic geometry than Fault I or Fault II. The seismic moment obtained using slip-grac-III is about 5.0×10²² Nm (equivalent to Mw=9.07), similar to that obtained using slip-grac-II and slip-grac-I. The slip distribution has maximum slip of about 22 m and a mean rake angle of 79°, consistent with other results from GRACE (Table 1). The slip pattern is distributed smoothly in the shallowest part of Fault III, but it is still clustered deeper than the pattern derived from GPS data (Figure 4a). Figure 8a shows a region where slip exceeds 10 m, extending along the strike of the fault for approximately 500 km. The region where slip exceeds 15 m extends for about 350 km, meaning that the slip area obtained using slip-grac-III is slightly wider than that obtained using cGPS data (which predicted a slip area of about 300 km along the strike (Zhou X et al., 2014). However, the slip area is more compact than the areas obtained using slip-grac-I and slip-grac-II. This can be seen in the high slip notches surrounding the epicenters. Unlike slips obtained by cGPS inversion, a large amount of slip smaller than 15 m was distributed at the south fault. Compared to the results of slip-grac-I and slip-grac-II, slip-grac-III results more closely surround the JMA epicenter. In addition, large ruptures occur around the hypocenter, a result in agreement with results inverted from GPS data. The centroid predicted using slip-grac-III, denoted by the triangle in Figure 7a, is closer to the JMA epicenter than the centroid predicted using slip-grac-II (Table 1). The centroids of the three models approximate those found by Cambiotti and Sabadini (2013). The centroid depth of 17.4 km predicted by slip-grac-III is slightly shallower than the centroid depths predicted by models I and II, is perfectly consistent with the centroids obtained by Han SC et al. (2011) and Cambiotti and Sabadini (2013), and is deeper than the centroid derived from GPS data. These results are consistent with the results of the sensitivity test described in Section 4.1 (Figure 3).

By using the slip distribution on Fault III, we predict the gravity changes and their residuals comparing with GRACE observations, as shown in Figure 8(b, c). Modeled gravity changes from -10.5 to 3.4 μGal account for ~72% of the GRACE data; this result is slightly improved compared to results obtained using slip-grac-I and slip-grac-II. Model results and observations show good agreement; the gravity change in the proximity of the dipolar pattern is particularly well reproduced, with a maximum residual of 1 μGal appearing in the region of negative gravity changes in Figure 8c.

We can use the length along the strike of slip over 15 m to quantify the spatial concentration of the slip inversion. The lengths of slip over 15 m for the slip-grac-I, slip-grac-II and slip-grac-III are 375 km, 425 km and 275 km, respectively. Similarly, the slip-cGPS-I, slip-cGPS-II, and slip-cGPS-III have spatial concentration distances of 300 km, 300 km, and 250 km, respectively. Although GRACE and cGPS behave similarly, portraying the tendency to concentrate slip for realistic strike changes in agreement with Slab1.0, only slip-grac-III matches slip-cGPS-III in terms of degree of slip concentration around the epicenter. The spatial concentrations, 275 km and 250 km, respectively, differ by only 10%, rather than by 25% or 42%, as was found for the less realistic fault models.

Figure  8.  Similar to Figure 5, slip-grac-III (a) inverted from co-seismic gravity changes estimated from GRACE data (Figure 1c). Co-seismic gravity changes modeled by the inverted slip-grac-III and associated residuals are shown in (b) and (c).

DownLoad: Full-Size Img PowerPoint

4.3.4   Comparison of GRACE and GPS results

Figure 9 shows that displacements predicted using slip-grac-I and slip-grac-III agree better with observations for both horizontal and vertical components than those predicted by slip-grac-II because their seismic moments agree better with those obtained using cGPS data. In particular, we observe a general improvement for both horizontal and vertical displacements in the region of (~140°E–141.5°E, ~36°N–38°N) in terms of amplitude and azimuth. To see this difference, compare Figure 9(a, b, e and f) for GRACE slip models I and III, with panels (c) and (d).

We use the root-mean-square (RMS)

to statistically elucidate differences between the model results and observed results. The RMS errors are 39.2, 58.8, and 45.1 cm in the horizontal and 13.6, 21.1, and 14.3 cm in the vertical displacements for models slip-grac-I, slip-grac-II, and slip-grac-III, respectively, while they are 2.33, 2.28, and 2.31 cm in horizontal and 2.71, 2.66, and 2.57 cm in the vertical components for model slip-cGPS-I, slip-cGPS-II, and slip-cGPS-III, respectively. We expect results based on cGPS data to exhibit smaller displacement RMS values than results based on GRACE data since cGPS is designed to sample this kind of data, whereas GRACE is designed to sample mass readjustment. The fault slip models slip-grac-I and slip-grac-III perform better than slip-grac-II, but the RMS error is greater than those associated with cGPS and teleseismic waveforms (Zhou X et al., 2012). This indicates that, due to the limitation of the temporal-spatial resolution of GRACE, the ability of GRACE to detect crustal deformation associated with a megathrust event is inferior to our ability to detect such deformation using waveforms. It is important to note that, because of the more compact pattern of slip-grac-III, we obtain better fits for both horizontal and vertical slip in the Japan region facing the epicenter, compared to fits obtained for slip-grac-I. We conclude that slip-grac-III performs better than the flat fault model at the single cGPS site with the most significant displacement.

Figure  9.  Comparison of co-seismic (a, c, e) horizontal and (b, d, f) vertical displacements between GPS observations (black arrows) and slip model predictions (red arrows) resulting from the slip model inverted from GRACE data (portrayed by the thin contour lines in Figure 5a, Fig 7a, and Figure 8a, where each contour represents 5 m). Displacements, from top to bottom, are predicted by slip-grac-I, slip-grac-II and slip-grac-III, respectively. The thick solid line with black triangles indicates the trench.

DownLoad: Full-Size Img PowerPoint

4.4   Joint Inversion

As shown in previous sections, the slip distributions inverted from the sole GRACE co-seismic gravity change and those measured by GPS do not completely match. To explain both cGPS and GRACE data, we need to jointly invert for the slip, combining the two datasets. In the joint inversion, we must determine the weight ratio of both datasets because of the discrepancies between the slip models, which are ultimately caused by the different nature of the two datasets (associated with deep masses or surface displacements, respectively). If we give a large weight to GRACE, the inverted slip will model the gravity changes well, but match the GPS measurements poorly, and vice versa. In the joint inversion, we fix the cGPS data and tune only the weight of GRACE data.

The results shown in Figure 10 indicate that as the weight of the gravity change measured by GRACE increases, the slip model produces a more extensive distribution with higher peak slip. The higher the gravity change, the better the fit between the modeled gravity change and the GRACE observations. However, the modeled displacements are in a better agreement with GPS measurements as the weight of inland GPS data increases (Figure S2). Taking the misfits of both datasets into consideration, we conclude that the ideal weighting ratio between inland GPS and GRACE data is 1:10 (slip-cGPS&grac), which produces a maximum slip of about 43 m and a mean rake angle of 74° (Table 1). The seismic moment is about 4.8×10²² Nm (equivalent to Mw=9.05), which is smaller than those obtained using slip-grac-II and slip-grac-III. The model can explain ~65% of the GRACE gravity amplitude. Meanwhile the RMS errors are 5.3 cm in the horizontal displacements and 5.1 cm in the vertical displacements, as determined by comparing the modeled and the observed cGPS data. This illustrates that the slip model slip-cGPS&grac is more localized and agrees better with GPS observations than the GRACE gravity change.

Although the cGPS data can improve the resolution and shrink the slip area, it cannot distinguish the slip near the trench (Zhou X et al., 2014). To improve the resolution close to the trench, co-seismic displacements of 7 GPS-A stations are involved in the joint inversion (Sato et al., 2011; Kido et al., 2011). The prediction of concentrated and shallow slip patterns close to the epicenter is most accurate for slip inversion over Fault III. Therefore, we invert for the slip distribution only within the frame of Fault III. Using the same technique, we obtained a joint slip model from cGPS, GPS-A, and GRACE data (Figure 11). The relative weights of the cGPS, GPS-A, and GRACE data, 1, 10, and 10 respectively, were selected following many tests. This optimal joint-slip model is denoted as slip-joint. The seismic moment of the inverted model is 4.8×10²² Nm, equivalent to M_W=9.05, similar to that obtained using slip-cGPS&grac. Compared with joint inversion results obtained using inland GPS and GRACE data, the use of GPS-A measurements improves the accuracy of the slip near the trench. Figure S3 shows that the displacements modeled using slip-joint are well in agreement with GPS observations. The joint-slip model can explain ~60% of the gravity change measured by GRACE—an improvement on results obtained using GPS data alone.

Figure  10.  The co-seismic slip models, as jointly inverted from GRACE gravity changes and inland GPS displacements. From top to the bottom, the left column shows slips obtained using different relative weights on GRACE data during inversion (5, 10, and 20). The middle column shows the GRACE gravity changes modeled by the three inverted slip models. The right column shows the residuals. We prefer the result obtained using a GRACE data weight of 10 (the middle) as the optimal joint slip model (denoted by slip-cGPS&grac). The centroid of the slip model slip-cGPS&grac is marked by a blue triangle.

DownLoad: Full-Size Img PowerPoint

Figure  11.  Slip model slip-joint (a) inverted from co-seismic gravity changes estimated from cGPS, GPS-A, and GRACE data. The co-seismic gravity change modeled by (a) and associated residuals are shown in (b) and (c), respectively.

DownLoad: Full-Size Img PowerPoint

5.   Discussion

In addition to GPS, InSAR, and seismic waveform, the temporal gravity field measured by the GRACE satellite mission can be used to retrieve the heterogeneous co-seismic slip distributions of megathrust earthquakes, regardless of GRACE’s lower resolution of about 300 km.

Due to the different temporal and spatial resolutions of GPS and GRACE data, it is reasonable to expect some discrepancies between the inverted slip patterns. Our results show that the gravity change is always underestimated by slip models based on cGPS co-seismic displacements, which are measured over the minutes before and after the event. Quantitatively, the slip models based on cGPS data for faults I, II, and III explain at most only 50% of the gravity change observed by GRACE. Similarly, although more seriously, none of the slip model results based on GRACE data alone (slip-grac-I, slip-grac-II, and slip-grac-III) fully match the inland GPS data (Figure 9).

In particular, slip-grac-I and slip-grac-II overestimate both horizontal and vertical displacements on the Japan islands, especially in the southeastern region of Honshu (~139.5°E–141°E, ~35°N–38°N). This is likely due to the GRACE inverted slip distribution extending to the southern fault with respect to the GPS inversion (Ozawa et al., 2011; Loveless and Meade, 2011; Zhou et al., 2014). Moreover, the azimuths of the horizontal displacements in northern Honshu are severely affected by the slip of the north fault, generally deviating to the north compared to observations. Slip-grac-III performs better at cGPS sites facing the epicenter, due to the arcuate subduction arc of Slab1.0.

After the main shock, there were three strong aftershocks larger than M_W 7.0 in one month. Although these events cause significant deformation, modeling using the Global Centroid Moment Tensor (GCMT) shows that the effects are an order of magnitude smaller than those caused by the main shock (Figure 12). The discrepancies of the slip distributions between GRACE and GPS inversions should not result from the large aftershocks. The detailed reasons for this will be discussed in the future.

The two joint inverted models slip-cGPS&grac and slip-joint not only estimate the inland GPS measurements properly but also match the gravity changes better than inverted models based on GPS data alone. An average rake of 74° for both models, coincident with Dai CL et al. (2014), has a small deviation from the pure dip-slip. The slip-joint model shows that the loci of GPS-A observations have upward slips. Taking the single GRACE inversion into account, we can infer that the positive gravity changes observed by GRACE require slip on the southern part of the fault as well as a strike-slip component. As is the case for point source inversion from GRACE (Han SC et al., 2011; Cambiotti and Sabadini, 2012; 2013; Dai CL et al., 2014), our work shows that the gravity change over the full space significantly increases the seismic moment of the event. It is worth noting that the slip-joint model exhibits a greater seismic moment than that resulting from the use of GPS data alone. Its predicted co-seismic displacements are coincident with measurements. Therefore, GRACE gravity change data, inverted from heterogeneous and realistic fault slip distributions as shown here, can be used to more accurately determine the seismic moment of a megathrust event, even given a poor distribution of geodetic stations.

Relative to the single cGPS inversion, the slip model jointly inferred from cGPS and GRACE data results in a slightly higher seismic moment and a larger maximum slip value. The maximum slip resulting from cGPS data is about 30 m, whereas it is about 45 m for the slip model involving the combination of two datasets. This is due to the contribution of GRACE data, which contains information on the gravity field response to the earthquake, especially the mass migration at the ocean. Although the maximum slip inverted from GRACE data alone is about 20 m, much smaller than that inverted using cGPS data, the larger slip area means that it still has the higher seismic moment. Thus, it is reasonable to state that as the weight given to GRACE increases, the maximum also increases.

Due to the strip effect inherent to the GRACE observation system, we must smooth the co-seismic gravity changes using a spatial filter technique. Although this improves the signal to noise ratio of the long-wave length gravity field, it also impacts the fault slip inversion, especially affecting the spatial resolution. To obtain a slip distribution from GRACE comparable to that obtained from GPS data, a higher spherical harmonic degree and lighter filtering will be necessary. Combining GRACE and GOCE without any filter may be an effective way to address this issue (Fuchs et al., 2016).

Some studies on continental intraplate earthquakes have shown that near-field geodetic measurements such those made by GPS and InSAR, as important as seismological data (Chen JW et al., 2017), are very sensitive to complex fault geometries. For example, Shen ZK et al. (2009) and Wang Q et al. (2011) used several fault segments with complex geometries including variable strike and dip and a mechanical decoupling to model the 2008 Wenchuan (M_W 7.9) earthquake. This allowed them to invert for the slip distribution from InSAR and GPS measurements. For the 1992 Landers M_W 7.3 earthquake, an event involving dextral strike ruptures, Freymueller et al. (1994) considered five simple faults with different strikes as determined on the basis of surface rupture to fit the GPS data. Similarly, for the 1999 Hector Mine (M_W 7.1) event—which occurred on a high-dip fault—Jónsson et al. (2002) considered nine segments with variable strike and dip values using simulated annealing methods based on the minimum fitting error between the model and existing InSAR and GPS data. GRACE data will be sensitive to fault strike and dip angles when a point source inversion is carried out (e.g., Wang Q et al., 2011; Han SC et al., 2011; Cambiotti and Sabadini, 2012, 2013). However, our findings indicate that GRACE data are not very sensitive to the fault geometry; this is probably because the inverted slip model over an approximate fault resembles the inverted slip model over a realistic fault. There is no contradiction between the two findings because the relatively wide initial search range was adopted in the point source inversion while the slight variations among Faults I–III are investigated in this paper. To some extent, our joint slip models approach the accuracy of those derived using GPS and seismic waveforms.

Figure  12.  The gravity changes caused by the three aftershocks. Two occurred on March 11, 2011 (M_W 7.9 and M_W 7.6) and one occurred on April 7, 2011 (M_W 7.1). These are modeled by dislocation theory using seismic sources from GCMT.

DownLoad: Full-Size Img PowerPoint

6.   Conclusion

Previous GRACE studies on the seismic source parameters of the 2011 Tohoku-Oki earthquake (centroid, rake angle, strike, dip, and seismic moment) are based on point source or homogeneous slip fault schemes. In this study, we make the model more realistic by considering a finite fault and heterogeneous slip inversion based on spatial sampling of GRACE gravity data. We focus on the characteristics of slip inversion from GRACE gravity data and its role in joint inversion using GPS displacements.

Due to the ~300 km spatial resolution of GRACE, the inverted slip distribution is smoother than it is for other datasets (such as those derived from GPS, seismic waveform and tsunami data) and is not severely affected by the dip of the fault. If an appropriate weight ratio based on the data resolution is selected, the slip inverted from the combination of cGPS and GPS-A and GRACE data can closely reproduce the GPS displacements. GPS data not only improve the spatial resolution but also make the slip area more compact relative to the GRACE. When we include GRACE gravity changes in the model, the seismic moment of the event increases from ~3.5 to 4.8×10²² Nm. This is due to the longer sampling frame of GRACE relative to GPS.

Table  1.  Comparison of seismic source parameters for the 2011 Tohoku earthquake

Model	Data	Centroid			Rake (°)	Strike (°)	Dip (°)	Seismic moment (×1022 Nm)
		Longitude (°)	Latitude (°)	Depth (km)
Han SC et al.(2011)	GRACE L1B KBR	―	―	17	82	196a	10.5	5.4
Wang L et al. (2012c)	GRACE CSR RL04	―	―	―	88a	203a	10a	4.59
Cambiotti and Sabadini(2012)	GRACE GFZ RL04	142.78	37.83	17	87	―	―	5.3
Cambiotti and Sabadini(2013)	GRACE GFZ RL04	143.47±0.46	37.75±0.46	16.1	89.7±8.5	201.9±7.8	12.2±3.1	6.32±2.94
Dai CL et al.(2014)	GRACE CSR RL05	142.17	37.53	―	77.0	201a	10a	6.43
slip-cGPS-Ib	cGPS	142.84	38.10	15.0	94	201a	11a	3.55
slip-cGPS-IIb	cGPS	142.80	38.14	15.7	95	201a	variable	3.54
slip-cGPS-IIIb	cGPS	142.72	38.15	18.1	85	variable	variable	3.00
slip-grac-Ib	GRACE CSR RL05	142.46	37.70	20.0	80.9	201a	11a	4.94
slip-grac-IIb	GRACE CSR RL05	142.34	37.62	19.8	79.9	201a	variable	5.14
slip-grac-IIIb	GRACE CSR RL05	142.46	37.70	17.4	79.2	variable	variable	5.00
slip-cGPS&gracb	GRACE CSR RL05 & cGPS	142.78	37.88	16.7	74.0	variable	variable	4.79
slip-jointb	GRACE CSR RL05 & cGPS & GPS-A	142.90	37.91	14.5	74.0	variable	variable	4.81
Notes: a fixed parameter; b source parameters are obtained by weighted average over the finite rectangle fault, where slip on each subfault is used as the weighting parameter.

 | Show Table

DownLoad: CSV

Acknowledgments

We thank Dr. Shuhei Okubo of ERI, University of Tokyo, for his discussion on data processing. We are grateful to two anonymous reviewers for their helpful and constructive comments, which greatly improved our manuscript. The GPS data were provided by the ARIA team at JPL and Caltech. All original GEONET RINEX data were provided to Caltech by the Geospatial Information Authority (GSI) of Japan. This research was supported financially by the National Natural Science Foundation of China (No. 41574021, 41474059, 41331066, 41774088, and 41174063), a research grant from the Institute of Crustal Dynamics, China Earthquake Administration (No. ZDJ2017-23), the CAS/CAFEA International Partnership Program for Creative Research Teams (No. KZZD-EW-TZ-19), the Key Research Program of Frontier Sciences CAS (Chinese Academy of Sciences) (QYZDY-SSW-SYS003), the SKLGED foundation (SKLGED2014-1-1-E), the GOCE Italy Project (the Italian Space Agency and the European Space Agency Endorsement), and the China Postdoctoral Science Foundation (No. 133014).

Appendix A. GRACE Time-Series Analysis

In this study, to fit the time series of gravity changes we use a scheme similar to that of De Linage et al. (2009), but without an exponential term.

where the left side represents the observation data and the right side represents the model. ${a_0}$ is a constant; the second term is the long-term trend; the third term is the periodical signal over three overlapping periods—annual, seismic annual, and S2 tide wave (161 days). $b$ is the amplitude; $\omega$ is the angle frequency; $\varphi$ is the phase; $H$ is the Heaviside function, which represents the co-seismic jump in the fourth term; ${t_0}$ is the original time of earthquake; the fifth term represents the post-seismic decay; $\tau$ is the relaxation time. If the relaxation time is fixed, the linear Least Square can be used to estimate the coefficients in equation (A1) and fit the time series.

Appendix B. Observation Equation

Unlike in the Slepian function, we here use a spatial sampling to model gravity change. First, we compute gravity changes and vertical displacements in a high-resolution space over the study region. These are can be represented as

where $g_{(\theta ,\lambda )}^{ solid}$ and $h_{(\theta ,\lambda )}^{ solid}$ are the gravity change and vertical displacement of the solid Earth’s surface; ${{{G}}_g}$ and ${{{G}}_h}$ are corresponding Green’s functions; ${{s}}$ is the slip vector. There are many methods for computing theoretical Green’s functions of gravity change and vertical displacements, including the homogeneous half-space (Okada, 1985; Okubo, 1992), layered half-space (Wang et al., 2006) and spherically layered models (Sun WK and Okubo, 1993; Sun WK et al., 2009). Herein we adopt the spherically layered elastic dislocation theory (Sun WK and Okubo, 1993) to compute ${{{G}}_g}$ and ${{{G}}_h}$ on the PREM Earth model [Dziewonski and Anderson, 1981). It is worth noting that ${{(}}\theta {{, }}\lambda {{)}} \in {\Omega _0}$ , which is a high-resolution space gridded in 0.1° increments.

Ocean water will redistribute over the deformed crustal surface after an earthquake occurs, causing additional change to the gravity field. This effect has the same magnitude as that of changes in the solid Earth and cannot be ignored (De Linage et al., 2009). In order to consider a complete solution, including attraction and loading, it is necessary to solve a sea level equation (Broerse et al., 2011). However, the loading effect is a second-order effect that is significantly smaller than the direct attraction effect. This allows us to use an infinite plane model (Bouguer layer) to approximate gravity changes due to the ocean water redistribution:

where $g_{(\theta ,\lambda )}^{ total}$ is the total gravity change; the second term of the right is the Bouguer correction (Heiskanen and Moritz, 1967); $G$ is the universal gravitational constant; ${\rho _{ o}}$ is the ocean water density, and ${O_{(\theta ,\lambda )}}$ is the ocean function:

Next, we expand the total gravity change on the surface using spherical harmonics (Heiskanen and Moritz, 1967) as follows

where $\bar C_{lm}^{{{model}}}$ and $\bar S_{lm}^{{{model}}}$ are Stokes coefficients, ${\bar P_{lm}}(\cos \theta )$ is normalized Legendre function, and $l$ and $m$ are the degree and order of spherical harmonics, respectively. We then obtain the gravity change on GRACE space ${\Omega _G}$ gridded by 0.5°× 0.5° by truncation (same degree as GRACE) and summation. The observation equation can be written as

where $g_{(\theta ',\lambda ')}^{ GRA}$ is the GRACE-observed gravity, $(\theta ',\lambda ') \in {\Omega _G}$ , $W_{lm}^{ C}$ and $W_{lm}^{ S}$ are coefficients of the DDK filter (Kusche, 2007), δ ${\Omega _k}$ is the area of a single grid cell, $N\!\!f$ is the number of subfaults, $N\!d$ is the maximum degree (e.g. 60 for CSR, 90 for GFZ), $N$ is the number of grids of high resolution space ${\Omega _0}$ , $e$ is the difference between each observation and the model. If using an isotropic Gaussian filter instead, equation (B5) can be reduced to

where ${W_l}$ is the coefficient of the Gaussian filter, and ${\psi _k}$ is the angular distance between $({\theta _k},{\lambda _k})$ and $(\theta ',\lambda ')$ . The addition theorem is used for equation (B6) (Heiskanen and Moritz, 1967).