Geodetic data inversion to estimate a strain-rate field by introducing sparse modeling

We observe active seismicity and crustal deformation in subduction zones. Since earthquake occurrences are closely related to strain accumulation, it is important to accurately estimate a strain-rate field. Many studies have estimated spatially continuous strain-rate fields from spatially sporadic geodetic data such as GNSS (Global Navigation Satellite System). However, localized strain rates near fault zones have tended to be underestimated, because most studies have applied a smoothness constraint (e.g., Okazaki et al., 2021, EPS). To overcome this difficulty, we introduce sparse modeling into the estimation of a strain-rate field. In this study, for simplicity, we consider the anti-plane strain problem.

We firstly express a velocity field by the superposition of cubic B-spline functions. Then, considering that a strain-rate field is smooth in most areas but can change abruptly in a narrow zone such as a fault zone, we impose both the sparsity constraint and the smoothness constraint of strain rates, which are expressed by the L1-norm and the L2-norm of the second derivatives of the velocity field, respectively. The relative weights of these terms are specified by two hyperparameters; the optimal values of which are determined by using the leave-one-out cross-validation method. We obtain the optimal values of the expansion coefficients of the cubic B-spline functions by minimizing the objective function, which consists of the terms of data fitting, the sparsity constraint, and the smoothness constraint.

To investigate the validity and limitation of the proposed method, we conduct synthetic tests, in which we consider an anti-plane strain problem due to a steady slip on a buried strike-slip fault. As a result, we find: (1) regardless of the locking depth of the fault, the proposed method reproduces localized strain rates near the fault with almost equal or better accuracy than the L2 regularization method, which imposes only the smoothness constraint, (2) the advantage of the proposed method over the L2 regularization method is clearer when fewer observation points are available, and (3) the proposed method can be applied when observation errors are small.

Next, we apply the proposed method to the GNSS data across the Arima-Takatsuki fault zone, which is one of the most active strike-slip faults in Japan. The proposed method estimates about 1.0×10^-8/yr faster strain rates near the fault zone than the L2 regularization method, which corresponds to a 20-30% greater strain-rate concentration. The faster and more concentrated strain rates result in the estimation of a shallower locking depth. Fitting the analytical solution to the estimated strain-rate profile, we obtain the optimal values of locking depth and steady slip rate as 11 km and 4 mm/yr for the proposed method, while 17 km and 5 mm/yr for the L2 regularization method. Since the former is closer to the depth of D90, 12-14 km (Omuralieva et al., 2012, Tectonophysics), above which 90% of earthquakes occur, this result suggests that the proposed method estimates a more realistic locking depth than the L2 regularization method.