Introduction of covariance components in slip inversion of land and seafloor geodetic data and application to slip deficit rate estimation in the Nankai Trough subduction zone

When spatial distribution of observation stations has bias in geodetic slip inversion, modeling errors in the inversion scheme may result in significant unnatural short-wave components in estimated slip distribution, which overfit to data. Combined use of both land and seafloor geodetic data in slip inversion often leads to such situations. To avoid overfitting, I proposed a method to incorporate proper covariance components in the data covariance matrix originated by modeling errors in the geodetic slip inversion. Because linearity of the inversion problem is retained, widely-known approaches to introduce prior constraints to the inversion, which assumes linear inversion, are easily applied to the proposed method. Synthetic tests showed that the introduction of covariance components allowed for the estimation of slip distributions closer to the true one, avoiding overfitting to geodetic data in biasedly distributed observation stations, compared to a conventional approach that does not introduce covariance components. I applied the proposed method to the estimation of the slip deficit rate in the Nankai Trough subduction zone, using geodetic data of displacement rates provided by land GNSS stations and seafloor GNSS-Acoustic stations. The proposed method estimated a reasonably smooth distribution of slip deficit rate compared to the conventional approach. Spatial distribution of residuals in the displacement rates suggested that the proposed method successfully avoids overfitting to the block motion at the south of the Median Tectonic Line, which the physical model I assumed here cannot describe.


Introduction
The development of the observation technique introduced seafloor geodesy such as GNSS (Global Navigation Satellite System)-Acoustic (GNSS-A) techniques and ocean bottom pressure (OBP) gauges to geodetic slip inversion, which was previously performed mainly based on GNSS data. The introduction of seafloor geodesy has drastically increased the resolution of inversions and improved the constraints on coseismic slips (e.g. [13]), afterslips, post-earthquake deformation (e.g. [24,15,1]), slow slip events (e.g. [29,30]), slip deficit rates (e.g. [14,37]), and plate motion itself (e.g. [4]). However, the combined use of land and seafloor geodetic data often results in the spatially biased distribution of the stations, because seafloor observation is usually more sparse and is located in more limited regions than land observations due to the relatively high cost of installation and operation. If such bias in spatial distribution is not properly considered in the inversion scheme, overfitting to the data sometimes causes unnatural short-wave components in the estimated slip distribution. To avoid this, larger weights are usually put on seafloor data (e.g. [10,13,7]), or a downsampling of the land stations is carried out (e.g. [37]). However, it is generally difficult to objectively determine the parameter for weighting and downsampling.
In the studies of seismic waveform inversions, [32,33] pointed out that the introduction of proper covariance components in the data covariance matrix, which is often neglected in geodetic slip inversions, can avoid overfitting to biasedly-distributed data. They discussed the sampling rate of the waveform data; they estimated the covariance of data in the time series based on the propagation of modeling errors due to anelastic attenuation of seismic waves and succeeded in the coseismic slip inversion that does not depend on the sampling rate of the waveform data. The same idea should be applicable to the slip inversion of geodetic data with spatial bias by associating covariance components to modeling errors in the inversion scheme.
As shown later, modeling errors in geodetic inversion mainly consist of errors in Green's functions and slip distribution. [6,21] introduced covariance components that depend on model parameters to make estimations, focusing on the error in Green's functions due to factors such as the uncertainty of elastic parameters and dip angles. They formulated a nonlinear inversion based on a Bayesian sampling algorithm, while the geodetic slip inversion is commonly formulated as a linear problem.
Here we propose an alternative approach that estimates spatial covariance between each geodetic data based on the Green's function itself. Because this method does not require model parameters to calculate the data covariance matrix, widely accepted linear inversion approaches such as [31] can be applied in the same way used by [32] for seismic waveform inversions.
This paper is organized as follows: Section 2 describes the mathematical formulation. Section 3 presents synthetic tests and compares the proposed method with a conventional approach. In Section 4, I apply the proposed method to the inversion of geodetic data of interseismic displacement rates and show that overfitting to the data is avoided in the estimation of the slip deficit rate when the proposed approach is used. Concluding remarks are provided in Section 5.

Calculation of data covariance matrix
Let us discuss a simple linear inversion problem u = Gp, where u, G and p are a vector for observation data, an observation (Green's function) matrix to relate the observation and unknowns, and a vector for unknown parameters, respectively. I assume Gaussian distribution of all stochastic variables. In many studies of geodetic inversions, the misfit vector ε = u − Gp is considered to obey the Gaussian distribution, i.e., ε ∼ N (0, E), where E is a diagonal matrix. In such cases, covariance components in the data covariance matrix are totally ignored. However, based on more detailed analysis, ε can be decomposed as where ε o = u − G true (p true ) and ε m = G true (p true ) − Gp are error terms due to observation and modeling, respectively. How to set ε o is relatively well studied [9, 14, ], and the use of a diagonal matrix for the covariance matrix is also expected to be a fairly good approximation. Hereafter, let us suppose the proper probabilistic distribution function (PDF) for ε o ∼ N (0, σ 2 o E o ) is known. [33,6,21] assumed that errors of Green's function is the main source of the error term ε m . If so, it is obvious that ε m has stronger spatial correlation between each component. Simply put, [33,6,21] assumed ε m = G true (p true ) − Gp = δGp, where δG is error of Green's function. Such a manner requires estimating error of Green's function and performing nonlinear inversion because the data covariance matrix changes with the model parameters p. Instead, in this study, I assume where p ′ true = G + G true p true , and G + and G true is a general inverse matrix of G and a matrix that linearizes the nonlinear operator G true , respectively. Here, I use an approximation that GG + ≃ I, where I is an identity matrix. Let us assume a simple PDF for δp = p ′ true −p, such as δp ∼ N (0, σ 2 s E s ), where I use E s = I, for instance. Then, due to the propagation law of uncertainty, Because G is usually a non-sparse matrix, σ 2 s GE s G T has non-zero covariance components. Based on the additivity of the normal distribution, the likelihood function of u can be written as where Note that covariance components are introduced to the likelihood function only by calculating matrix-matrix products using G without any additional information on the inverse problem. Furthermore, because the linear property of the inverse problem is not changed, the proposed approach is easily applied to existing geodetic inversion methods such as [31].
As mentioned above, Equation 5 assumes GG + = I, which holds when G's number of columns is greater than the number of rows. Otherwise the equality is just an approximation. The effect of this approximation on the estimation results will be discussed in Section 4. In addition, E s = I implicitly assumes G + G true ≃ I. I should note that an excessive difference between G and G true may reduce the estimation accuracy of the prior PDF.

Applicability to the inversion formulation based on Bayes' theorem and Akaike's Bayesian information criterion
I demonstrate the applicability of the proposed approach to an existing geodetic linear analysis method proposed by [31]. Their method incorporated prior constraints of the smoothness of slip distribution to the observation equations, constructing a Bayesian model with unknown hyperparameters. The values of the hyperparameters, which control the structure of the Bayesian model, are objectively determined from observation data by minimizing Akaike's Bayesian information criterion (ABIC) [2, ]. I can directly apply the likelihood function defined in Equation 7 to their approach. I first assume a prior PDF with respect to the smoothness of the slip distribution, where L is a matrix corresponding to a discretized Laplacian operator, P is the rank of L T L, and Λ is a P × P diagonal matrix that has non-zero eigenvalues of L T L as its diagonal components. Based on Bayes' theorem, the posterior PDF of the model parameters is written as where Optimal model parameters that minimize this functional is readily obtained once the hyperparameters are determined, as Here, I introduce ABIC, which is defined as where The necessary condition for ABIC to be minimized is which allows for ABIC being dependent on two hyperparameters (See Appendix A for detailed calculation of ABIC). I adopt the hyperparameters that minimize ABIC in the inversion, performing minimization using a grid search. Thus, I can determine the hyperparameters for the PDF for the proposed method by applying [31].

Synthetic test
In this section, I perform synthetic tests of geodetic inversion that include significant modeling errors and compare results obtained using the proposed method and a conventional one. When using the conventional method, I assume ε ∼ N (0, σ 2 E), where σ is a hyperparameter and E is a diagonal matrix. The approach of [31] is used elsewhere, determining the hyperparameter on the prior PDF of the parameter based on ABIC in the same way as described in Section 2.2. I show the advantage of introducing covariance components due to modeling errors in the geodetic inversion analyses.

Problem setting
I consider a dip slip fault in the elastic homogeneous half-space, which models a plate boundary. The fault geometry is shown in Figure 1 (a): the fault width W = 680 km, the fault length L = 340 km, and the dip angle θ = 40 • , respectively. I consider a synthetic test problem of estimating slip deficit distribution (SDR), which is targeted in the application of actual data in Section 4. SDR and geodetic data for inversion have the dimension of velocity. Otherwise, it is the same as in ordinary slip inversions. The true SDR distribution shown in Figure 1 (b) is inputted in Plane 1 and elastic responses in virtual observation points are used as the synthetic data. The virtual observation points are biasedly distributed, which mimics geodetic inversion using both land and seafloor stations. Here, I consider two patterns of observation distribution, the main difference of which is the location of the seafloor stations, as shown by red and green points in Figure 1 (b). I assume that 16 relatively isolated points located direcly above the main source area are seafloor stations and 800 other points land stations. Following the configuration of the real world problem presented in Section 4, I use the x, y and z components of the displacement rates in the land stations, and the x and y for the seafloor stations. I add artificial Gaussian noise to the computed displacements, for which the standard deviations are 2.5×10 −4 m/yr for the horizontal component and 5.0×10 −4 m/yr for the vertical component. The covariance matrix for the observation error, E o for the proposed method and E for the conventional method, is a diagonal matrix, for which the components are set according to these noise levels. To introduce modeling errors in the synthetic tests, I calculate Green's functions in Plane 2, a fault plane that has a perturbation dip angle dθ from the true one ( Figure 1 (a)). The size of subfaults h to input unit slip for calculating Green's functions is 20 km. I consider dθ = 2 • , 5 • and 10 • for comparison. I calculate both the synthetic data and Green's functions using the computation program of [20]. Using the synthetic data and the Green's functions, I perform inversion analyses using the proposed method and the conventional method and compare the estimated SDR distribution.

Results
When dθ = 2 • , the estimation results look almost identical, regardless of which method is used ( Figure  2). It implies that when modeling errors are small, neglecting covariance components does not cause significant differences in the estimation results. When dθ = 5 • , modeling errors become larger than in the case of dθ = 2 • . As a result, we observe significant differences in the estimation results ( Figure 3). The proposed method estimates a smoother SDR distribution, which is closer to the true one than the conventional method; the root mean square error (RMSE) of the estimated slip deficit distribution is 3.4×10 −3 m/yr for the proposed method and 7.0×10 −3 m/yr for the conventional one. Such a tendency becomes more obvious when dθ = 10 • ; many significant short-wave components, for which the spatial distribution seems strongly affected by the locations of observation stations, are seen in the estimation result of the conventional method, which are not included in that of the proposed one ( Figure 4). Although the larger error in Green's functions also degrades the estimation results of the proposed method with SDR peaks imposed to the region along the trench axis, the basic feature of the true SDR, which consists of two peaks, remains in the Figure 4 (a). The RSME of the SDR distribution in the case of the proposed method is 1.2×10 −2 m/yr, which is significantly smaller than that for the conventional method, 1.7×10 −2 m/yr. Such tendency strongly appears when a more biased distribution of observation stations is adopted for dθ = 10 • ( Figure 5). The amplitude of short-wave components in the SDR estimated by the conventional method is even larger in this case, whereas the proposed method still estimates a smooth distribution regardless of the location of the stations. A relatively large discrepancy between observed and simulated displacements shown in Figure 5 (c)(e) implies that overfitting to the data is avoided in the case of the proposed method.
In summary, the proposed method gives smoother solutions that have smaller RMSE to the true SDR, compared to the conventional method. When the modeling errors are small, the proposed method provides almost identical estimation results to that of the conventional method. These findings suggests that the proposed method succeeds in estimating the modeling error and PDFs, which produces favorable estimation results.

Application to slip deficit rate estimation in Nankai Trough subduction zone
The Nankai Trough subduction zone, located off the coast of southwest Japan, is well known for great earthquakes with recurrence intervals of 100-200 years [3, ]. To estimate the accumulated stress on the focal region in the present day, it is important to estimate the SDR on the plate boundary. The interseismic displacement rates in southwest Japan have been compiled by continuous GNSS observation provided by GNSS Earth Observation Network System (GEONET), which is managed by the Geospatial Information Authority of Japan [22, ]. In addition, seafloor geodetic observation using GNSS-A has provided interseismic displacement rates in the region directly above the megathrust [37, ]. While combined use of the land and seafloor geodetic data is promising in increasing spatial resolution of SDR estimation compared to the case solely using the land data, spatial bias in the distribution of observation stations may lead to a severe problem of overfitting as discussed in Section 1. To address this problem, I apply the proposed method to the estimation of SDR in the Nankai Trough subduction zone and discuss its validity.

Data and modeling setting
I use GNSS and GNSS-A data of the displacement rates, which were estimated by [18] with a common reference frame fixed to the Eurasian plate. The GNSS displacement rates were estimated for the east, As a result, 316 GNSS stations and 19 GNSS-A stations are used. I use homogeneous elastic halfspace to calculate Green's functions. It means that the displacement rate data is expected to include at least two major sources of modeling errors; one is displacement due to rigid block motions. The Median Tectonic Line (MTL) is a famous arc-parallel fault system in this region, in which the rightlateral strike-slip movement is related to the oblique subduction of the Philippine Sea plate. Past studies estimated that this movement imposes around 5 mm/yr motion in the southern block of MTL [17, 25, ]. Another factor is the 3D heterogeneous elastic structure that cannot be incorporated in the homogeneous elastic half-space. Plate boundary geometry is extracted from "Japan integrated velocity structure model version 1" (JVSM) [16, ], a unified velocity structure model for the Japanese Islands. Plate boundary is discretized by 4km-sized triangles. I use 3rd-order bicubic B-splines for basis function to discretize the slip deficit distribution following [31], central points of which are in 20 km intervals ( Figure 6). I locate 432 central points, so that they discretize the area of the plate boundary between 131.5 • E-139 • N, 31.5 • N-35 • N, and shallower than 60 km depth ( Figure 6). Because we estimate SDR in two components (the direction of plate convergence and the one perpendicular to it), the number of unknowns is 864. Response function due to a triangular fault in the elastic homogeneous half-space is calculated based on an analytical expression of elastic response due to an angular dislocation proposed by [5].
In the inversion scheme, I assume E o = I for the proposed method and E = I for the conventional method for simplicity.

Results and Discussions
Figure 7 (a) shows the SDR distribution estimated using the proposed method. A smooth solution, which is not strongly correlated with the positions of geodetic stations, is obtained. In contrast, SDR distribution estimated by the conventional method has significant short-wave components directly below the area in which geodetic stations are absent (Figure 7 (b)). The maximum slip deficit rate in this subduction zone is expected to be no more than the reported plate convergence rate, i.e., at most around 6 cm/yr [11, ]. This knowledge in plate tectonics is more consistent to the estimation result of the proposed method. It follows that the proposed method avoids overfitting to the factors that are not modeled in the scheme, which otherwise become sources of unnatural short-wave components in the estimated slip. ABIC is -9.658×10 3 for the proposed method and -9.541×10 3 for the conventional method, which suggests that the SDR model estimated using the proposed method is more preferred in terms of model selection as well. Here are the detailed discussions: Figure 8 compares the residuals of surface horizontal displacement rate between the proposed and conventional methods. Systematic patterns are observed in the result of the proposed method. Part of these residuals are expected to originate from the displacement rates due to factors that are not modeled in the scheme. In particular, we observe a significant change of residuals around MTL and the Beppu-Shibabara graben (BSG), west of MTL; significant displacement rates towards the west are seen in the land station at the south of MTL and BSG, but are not seen at the north. These residuals in the west direction in the southern block of MTL are at around the rate of 5 mm/yr, which is consistent to the implication on block motion in the southern block of MTL from past studies [17, 25, ]. This strongly indicates that our method succeeds in avoiding overfitting to the displacement rates produced by relative motion along MTL and BSG. Such systematic residuals around MTL and BSG are not seen in the results of the conventional one. However, we also see strong systematic patterns in other regions such as the Tokai area, for which I cannot suggest a reasonable explanation. While the proposed method is expected to reduce the effect of the modeling errors, it is, needless to say, also important to make efforts to reduce modeling errors such as the consideration of block motion [18, 19, ] and incorporation of 3D elastic heterogeneity and viscoelasticity [27, 19, ]. Figure 9 shows the resolution of the SDR estimations. The resolution I discuss here is the diagonal component of the model resolution matrix. In our formulation, the model resolution matrix R m is computed as:

Resolution of estimated SDR
The L-2 norm of the resolution for two direction of slips in each subfault is plotted. Interestingly, the resolution of the SDR distribution estimated by the proposed method is significantly low. This is reasonable because the larger weight put on the smoothness constraint, the lower resolution obtained in the estimation result, shown in Equation 16. This indicates that we cannot succeed in both avoiding overfitting and achieving high-resolution estimation when the problem includes significant modeling errors.

Data downsampling and data covariance components
One of the important issues in the combined use of land and seafloor data for conventional geodetic slip inversion is how to determine the weight on each kind of data; some previous studies place larger weight on seafloor data [13, ] or downsample the land stations [37, ] to seek a better balance between the land and seafloor data. One advantage of the proposed method is that introducing proper data covariance prevents us from determining weights on the data. Figure 10 shows estimation results with 100 randomly sampled land stations out of the original set of 316 stations. The slip distributions estimated using the proposed method are not sensitive to the number of land stations used in the estimation, whereas it is greatly affected when using the conventional method. Table 1 lists RSME in each case of estimation. In addition, downsampling of the land stations does not affect the RMSE much regarding both all and the seafloor data for the proposed method, while it affects significantly for the conventional one. It follows that the introduction of the proper data covariance matrix avoids the problem of determining weights on each kind of observation data in the cost function of an inversion problem. When all the stations are used in the inversion, G's number of rows is 316×3+19×2=986, which is greater than the number of columns, 832. This means GG + = I, which Equation 5 in Session 2 assumes, just approximately holds. In the case of downsampling, G's number of rows is 100×3+19×2=338, which is smaller than the number of columns, with which GG + = I holds. In spite of such a difference, the estimated SDR distributions in the two cases is similar, which suggests that the approximation GG + ≃ I in the former case does not affect much the estimation result.

Conclusion
I proposed a method to incorporate proper covariance components originated by modeling errors in geodetic slip inversion to avoid overfitting to geodetic data in biasedly distributed observation stations. Because this method retains linearlity of the inverse problem, it is applicable to widelyknown approaches to introduce prior constraint to the inversion, such as the method of [31]. I carried out synthetic tests of SDR estimation in a dip slip fault, in which I introduce modeling error in estimation using different dip angles to the true fault geometry. The proposed method provided smoother solutions that have smaller RMSE to the true SDR assumed in the synthetic test, compared to the conventional method. I applied the proposed method to estimation of SDR in the Nankai Trough subduction zone, using geodetic data of displacement rates provided by land GNSS stations and seafloor GNSS-Acoustic stations. The proposed method estimated a reasonably smooth estimation result compared to a conventional approach that does not introduce covariance components. Spatial distribution of residual in displacement rates suggested that the proposed method successfully avoids overfitting to block motion at the south of Median Tectonic Line, which the model that I assumed  here cannot describe. It was also shown that it is difficult to succeed in both avoiding overfitting and achieving high-resolution estimation when the problem includes significant modeling error. Estimation results using downsampled land data shows that the estimation results obtained using the proposed method are robust to a change in sparseness of the stations, which suggests that the proposed method does not require subjective determination of weight in each data. As indicated by the examples in seismic waveform inversions [32, 33, ], the proposed approach should be also effective in time-dependent geodetic slip inversion, such as post-earthquake slip inversion considering viscoelasticity [34, ]. Introduction of proper covariance components in time is also expected to allow for robust estimation results to irregular time intervals in time series of seafloor geodetic data from GNSS-Acoustic stations (e.g. [28,26]). However, increase of the number of unknowns in time-dependent inversion may lead to difficulty in calculating inverse of the covariance matrix, which requires further consideration in computation algorithm.
Although the proposed method can reduce the effect of modeling errors in estimation results, the best way to improve the inversion results is of course to reduce the modeling errors itself by revising the numerical model. Recent development of numerical simulation and high-performance computing methods is a key factor for this purpose [12, 8, ]. In the future work, more direct consideration of uncertainty in Green's function is desirable, which should be introduced by executing multiple computations of crustal deformation [35, 36, ].