Potency density tensor inversion of complex body waveforms with time- adaptive smoothing constraint

1Graduate School of Science and Technology, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan. 2Faculty of Life and Environmental Sciences, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan. 3Mountain Science Center, University of Tsukuba, Tsukuba, Ibaraki 305-8572, Japan. 4COMET, School of Earth and Environment, University of Leeds LS2 9JT, UK. 5Japan Agency for Marine-Earth Science and Technology, 3173-25 Showa-machi, Kanazawa-ku, Yokohama 236-0001, Japan. 6Disaster Prevention Research Institute, Kyoto University, Uji, Kyoto 611-0011, Japan. Corresponding author: Shinji Yamashita (ijnihsatihsamay@gmail.com)


Introduction
Large earthquakes are often accompanied by complex fault geometries and slip distributions. Elucidating these complexities is important not only for understanding the source process itself, but also for evaluating the stress field (e.g. Dreger et al., 2004;Nakamura et al., 2010) and future seismicity around the source region (e.g. Hayes et al., 2010;Meng et al., 2012). With the recent development of source process analysis methods, it has been proven that observed body waveforms are sometimes difficult to be explained by the seismic source model using simplified fault plane assumptions (e.g. Shimizu et al., 2020;Yamashita et al., 2021).
It has been challenging to estimate a complex rupture process of an earthquake.
There are two major methods to analyze the seismic source process using seismic waveforms: finite-fault inversion (FFI) and back-projection (BP). While the FFI methods can obtain the spatiotemporal distribution of slip rates on an assumed fault plane, it has a weakness in that the obtained result varies depending on the assumption of fault model and underground structure (Yagi & Fukahata 2011, Minson et al. 2013, Duputel et al. 2014, Ragon et al. 2018. While the BP methods can estimate the seismic radiation by stacking the observed waveforms corrected by the travel-time shifts and the polarity of P-waves, the results are affected by reflected phases , Okuwaki et al. 2015, Lay et al. 2018, and the theoretical resolution is lower than that of waveform inversion (Fukahata et al. 2014), making it difficult in principle to determine whether the signals included in the results are ghosts or not. Both the FFI and the BP methods easily distort the solution when applied to earthquakes with complex fault geometries. In the case of the FFI methods, the discrepancy between the true and assumed fault geometries causes a modelling error that biases the obtained solution (Ragon et al. 2018 and can lead to incorrect estimation of the rupture propagation direction ). In the case of the BP methods, the variation of the fault geometry sometimes leads to temporal changes of the polarity flips, making it difficult to This is a non-peer reviewed EarthArXiv preprint 5 obtain coherent signals. There is an urgent need to develop an analysis tool that can provide robust source processes even for earthquakes with complex fault geometries.
Recently, a novel FFI method inverting teleseismic P waveforms has been developed to estimate spatiotemporal rupture evolution with the variation of focal mechanisms, projected onto an assumed model plane . Because this FFI method expresses estimated slip amounts with focal mechanisms as potency density tensors (e.g. Ampuero, 2005) on the assumed model plane, we refer it to a potency density tensor inversion. This inversion method expresses a slip direction on the assumed model plane as a focal mechanism represented by a linear combination of the five basis doublecouple components (Kikuchi & Kanamori 1991) (Fig. 1), so resulting slip direction is not confined to the assumed strike and dip angles of the model plane. In addition, the Green's function of the teleseismic waveform is insensitive to small deviations in the assumed source location, but sensitive to the assumed fault orientation . Thus, the assumed model plane no longer requires to be compatible with the true fault geometry.
This method significantly reduces modelling errors originated from incorrect fault orientations, leading to robust seismic source modelling (Hicks et al. 2020. Since the Green's function error is explicitly introduced into the data covariance matrix following Yagi & Fukahata (2011), the slip direction is stably obtained without a non-negative slip constraint (Yagi & Fukahata 2011, which has been widely used in the conventional FFI methods (e.g. Hartzell & Heaton 1983, Das & Kostrov 1990).
The potency density tensor inversion method employs a smoothing constraint as prior information, such that the spatial Laplacian and the second-order time derivative of the potency-rate density function are close to zero . Hyperparameters that control the strength of the smoothing constraint are objectively determined by Akaike's Bayesian Information Criterion from observed data (ABIC; Akaike, 1980;Yabuki & Matsu'ura, 1992). However, since the same hyper-parameters are used for all the basis double-couple components , relatively This is a non-peer reviewed EarthArXiv preprint 6 strong smoothing is imposed on the component with larger amplitude , which makes the inversion solution overly smooth. To mitigate inappropriate strength of the smoothing constraint among the five basis double-couple components, Yamashita et al. (2021) developed formulation that gives a weight of each basis doublecouple component in proportion to its moment tensor solution (e.g. Global Centroid Moment Tensor solution; Dziewonski et al., 1981;Ekström et al., 2012). Although their smoothing constraint has been successfully applied to the 2018 Gulf of Alaska ) and 2020 Caribbean (Tadapansawut et al. 2021) earthquakes, it refers to a single moment tensor solution integrated with respect to time and is therefore inflexible for earthquakes where the focal mechanism changes drastically during the rupture, such as the 2002 Denali fault earthquake (Eberhart-Phillips et al. 2003). In addition, their smoothing constraint introduces the constant weight for each basis doublecouple component over the entire source time, although the amplitude of the potency-rate density function of each basis double-couple component is generally time-varying. As a result, when the amplitude was small, relatively weak smoothing is applied, and when the amplitude was large, relatively strong smoothing is applied, resulting in inappropriate smoothness of the solution in the time direction. In order to mitigate this problem, the weights that control the strength of the smoothing need to be time-adaptive as the amplitude changes.
Here, we propose a time-adaptive smoothing constraint that dynamically controls the smoothing hyper-parameters for each basis double-couple component.
Following the mathematical reviews of an observation equation for the potency density tensor inversion and formulations of the conventional ) and our proposed time-adaptive smoothing constraints, we perform a numerical experiment for a synthetic model with significant changes of the focal mechanism and moment release during rupture, and compare the conventional and time-adaptive smoothing constraints. The numerical experiment shows that although the overall focal mechanism variation of the synthetic model can be obtained with the conventional This is a non-peer reviewed EarthArXiv preprint 7 approaches, the time-adaptive smoothing constraint enables us to better reproduce the amplitudes of each basis double-coupled component, which mitigates over-and undersmoothing of slip and fault geometry changes in space and time. We then further evaluate the feasibility and effectiveness of the time-adaptive smoothing constraint on real earthquakes with complex fault rupture. For the applications, we chose the 2002 Denali fault, Alaska and the 2008 Wenchuan, China earthquakes. We finally discuss the advantages and limitations of our new inversion method with the time-adaptive smoothing constraint.

Method
The potency density tensor inversion method ) expresses the slip direction on the assumed model plane as a superposition of five basis double-couple components (Kikuchi & Kanamori 1991) (Fig. 1),introducing an error term of the Green's function explicitly (Yagi & Fukahata, 2011). An observed teleseismic waveform ! at a station in the potency density tensor inversion can be expressed as where "! is the calculated theoretical Green's function of the -th basis double-couple component, "! is the model error on "! (Yagi & Fukahata 2011), ̇" is the th potency-rate density function, '! is a background and instrumental noise, represents a position on , and * denotes the convolution operator in the time domain.
As a prior constraint, Shimizu et al. (2020) adopted the smoothing constraint for the potency-rate density function in space and time, respectively, which can be expressed as where ) and * are noises. Eqs (2) and (3) can be rewritten in vector form as follows This is a non-peer reviewed EarthArXiv preprint  1980, Yabuki & Matsu'ura 1992, Fukahata et al. 2003 *" where is a scaling factor and " is the absolute total potency of the -th basis doublecouple component, which is independently derived from the moment tensor solution. The introduction of extremely small standard deviations imposes strong smoothing on the potency-rate density function, and leads to the amplitude fluctuations in parts where the smoothing is relatively weak, resulting in an unstable solution. To avoid this, the minimum value of " was adjusted by an given threshold. In the following, we refer to this smoothing constraint as the moment-tensor-weighted smoothing constraint.
Although the moment-tensor-weighted smoothing constraint ) has been successfully applied to several earthquakes (Tadapansawut et al. 2021, the method only refers to a single moment tensor solution and is therefore inflexible for earthquakes with significant changes in focal mechanism during rupture. In addition, since the value of " in eqs (6) and (7) is not time-varying over the entire source time, excessively weak smoothing is applied to the part of a potency-rate density function with particularly small amplitude, resulting in an unstable solution. To solve these disadvantages, in this study, we propose a time-adaptive smoothing constraint that time-adaptively changes the values of the standard deviations referring to the potency-rate function for -th component " , which is obtained by integrating the potency-rate density function in space. Then, the covariances ) L " ( ) ( and * L " ( ) ( for the time-adaptive smoothing constraint in space and time, respectively, can be expressed where N is a scaling factor. Since the covariances of the Gaussian noises in eqs (8) and (9) depend on the potency-rate function determined from the model parameters, the inversion problem to be solved becomes non-linear. The potency density tensor inversion framework originally includes a non-linear problem because the data covariance matrix, which introduces the error term of the Green's function, is a function of the model parameters. To solve this non-linear problem, Shimizu et al. (2020) gave initial model parameters to evaluate the data covariance matrix and iteratively solved it by improving the model parameters, following Yagi & Fukahata (2011). Therefore, in the same way, we can introduce the time-adaptive smoothing constraint with nonlinearity into the potency density tensor inversion framework. We introduce an input potency-rate function " +,-.* obtained from the initial given model parameters in the first iteration and from the last obtained model parameters in subsequent iterations into eqs (8) and (9). The potencyrate function, obtained by spatially integrating the potency-rate density function, can take zero or extremely small values. As explained above, the introduction of such extremely small values for the standard deviations results in an unstable solution. Thus, to prevent the input potency-rate function from having extremely small values, we defined the input potency-rate function " +,-.* for -th iteration as follows where ",,0& is the estimated potency-rate function at the ( − 1)-th iteration, 123 is the maximum value of absolute ",,0& for all , and is a lower threshold. Note that in the first iteration ( = 1), ",4 is calculated from the given initial model parameters. The scaling factor N in eqs (8) and (9) is set so that the minimum value of N ( " +,-.* ( ) ( becomes 1. Following Yagi & Fukahata (2011) and Shimizu et al. (2020), the hyperparameters including ) ( and * ( for each iteration are determined by minimizing ABIC, and the inversion iterations are repeated until the normalized L2 norm between the model parameters of the input and the result becomes acceptably small. We have confirmed that the final solution of the inversion iteration is independent of the setting of the initial model parameters through the following numerical experiment using an input model of next part and applications to real earthquakes shown in the discussion section.

Numerical experiment
In order to evaluate the performance of the time-adaptive smoothing constraint compared to the unweighted ) and moment-tensor-weighted ) smoothing constraints, we performed a numerical experiment. We inverted  (Table   S1) and were added random Gaussian noise. In addition to the noise of the Green's functions, random Gaussian noises are added to the input synthetic waveforms as the background noise '! . We generated the input synthetic waveforms at the same 45 teleseismic stations (Fig. S2a) as those used in the analysis of the 2002 Denali fault earthquake, which is described later.
Because the potency density tensor inversion method allows for any types of faulting mechanism on the assumed model plane, the inversion can be performed using a horizontal model plane, where the modelling error of fault orientation is replaced by that of fault location, which is less serious in the Green's functions of teleseismic waveforms ). Here, we set a horizontal rectangular model plane of 220 × 60 km in length and width at a depth of 25 km covering the input faults ( Fig. 2a) and discretized it into space knots of 10 × 10 km in length and width. The initial-rupture point was put at the same point as the input. The potency-rate density function at each space knot was represented as a combination of the bilinear B-spline functions with a time interval of 1.0 s, which had a total duration of 15 s. The maximum rupture front velocity, which defines the rupture starting time at each space knot, was set to 3.3 km/s. We performed inversions with the unweighted, moment-tensor-weighted and time-adaptive smoothing constraints, respectively.
For the moment-tensor-weighted smoothing constraint, following the numerical experiment of Yamashita et al. (2021), " in eqs (6) and (7) was calculated by referring to the input total moment tensor with a lower bound of 10% of max ( " ). As seen in the paragraph above, since the input total moment tensor is a pure strike-slip faulting (inset of Fig. 2a), the moment-tensor-weighted smoothing constraint imposed relatively large standard deviations on the M1 and M2 components corresponding to the strike-slip components ( Fig. 1).
In the time-adaptive smoothing constraint, we set the lower bound in Eq. (10), and the initial model parameters to obtain the input potency-rate function for the first inversion iteration. We set the lower bound as 0.2 in the following analyses.
Sensitivity of the solution to the value of was evaluated with = 0.2 ± 0.1, but the selection of did not affect the solution (Fig. S3). For these tests, we gave a uniform slip model with M1 slip component as the initial model parameters (Fig. S4a). The slip of each space knot is triggered by the presumed rupture front and decays gradually during the rupture duration. Because the time-adaptive smoothing constraint obtains the initial input potency-rate function from the given initial model parameters, the final solution may depend on the initial model parameter settings. To investigate the dependence, we inverted the numerical models with three different initial conditions: a uniform pure strike slip (M1), a uniform reverse slip (M5), and a striped slip (Fig. S4). The obtained results do not show a major difference ( Fig. S4), indicating the independence of the final solution to the initial model parameter settings. In the following, we use the uniform M1 slip as an initial model parameter. Figure 2 shows the potency density tensor distribution obtained by the three types of smoothing constraints. Figure 3 shows the input and resultant moment rate functions of the M1, M2 and M5 slip components included in the input model, obtained by spatially integrating the potency-rate density function multiplied by the rigidity.
Although the model geometries are different between the input fault and the model fault and the quantitative measure of model discrepancy may not be rigorously made in that context, the comparison of the total seismic moment can be useful indexes of discrepancy from the input.
In the unweighted smoothing constraint (Fig. 2b), the obtained spatial focal mechanism distribution well captures the input fault geometry. On the other hand, the obtained total seismic moment was underestimated by 15% of the input.   solution (Dziewonski et al. 1981, Ekström et al. 2012) is located about 130 km east of the epicenter and shows a predominantly strike-slip faulting (Fig. 4). Aftershock seismicity ( ≥ 3.5, 3-9 November 2002) (U.S. Geological Survey Earthquake Hazards Program 2017) is concentrated around the ruptured three faults (Fig. 4).
Given the distributions of surface rupture traces (Koehler et al. 2012) and aftershocks (U.S. Geological Survey Earthquake Hazards Program 2017) (Fig. 4), we set up a horizontal model plane of 336 × 96 km in length and width with a strike of 296°, which is the same as the strike angle of the north-dipping nodal plane of the GCMT solution (Dziewonski et al. 1981, Ekström et al. 2012. We adopted the epicentral location determined by the USGS NEIC for the initial-rupture point at 63.514°N and 147.453°W ( Fig. 4). Since the aftershock distribution is concentrated in shallow depth of 10 km, suggesting that the rupture has not reached deep underground (Eberhart-Phillips et al.

2003)
, we placed the model plane at a depth of 5 km. The potency-rate density function along the model plane was expanded using a bilinear B-spline with a spatial interval of 12 km and a temporal interval of 0.8 s with a total duration of 44 s. The potency-rate density function was estimated up to 100 s from the origin time (e.g. Ji et al. 2004, Asano 2005. The maximum rupture front velocity was set to 5 km/s to account for the possibility of supershear rupture propagation (Dunham & Archuleta 2004, Ellsworth et al. 2004, Frankel 2004, Asano 2005, Walker & Shearer 2009).
The obtained moment rate function has peaks at 9.6, 32.0 and 59.2 s (Fig. 5a) (Fig. 5b). The time evolution of the potency-rate density along the strike direction shows that the three peaks of the moment rate function correspond to large potency-rate density patches near the epicenter, around the intersection of the Denali fault and the Trans-Alaska Pipeline System (TAPS), and just before the junction between the Denali and Totschunda faults, respectively (Fig. 5c). These spatial distributions are consistent with the estimated location of the three sub-events identified in the strong motion studies (Eberhart-Phillips et al. 2003, Frankel 2004. The snapshot shows that the rupture initiated around the Susitna Glacier fault (Fig. 5d). The moment tensor solution obtained by integrating the potency-rate density tensors in time and space for each time window are shown in the inset of Fig. 5d. The initial rupture was dominated by reverse faulting and lasted until about 20 s from the origin time (Fig. 5d). Then, the rupture propagates unilaterally from west to east along the Denali-Totschunda fault with predominantly strike-slip faulting (Fig. 5d). The strike direction of the nodal plane derived from each moment tensor matches the orientation of the Denali-Totschunda fault from 20 to 80 s with a gradual rotation of 16° clockwise from 280° to 116° (Fig. 5d), which is consistent with the negative peak of the M2 component lagging behind the positive peak of the M1 component in the moment rate function (Fig. 5b). The synthesized waveforms obtained our source model explain the observed waveforms well (Figs. 5e and S5b).
The centroid migration velocity of the potency-rate density decelerates around the Denali fault-TAPS intersection but accelerates before and after it (Fig. 5c). The estimated average centroid location migration speeds are 3.16 km/s for 10-40 s, where the centroid is located west of the TAPS. Then, the migration velocity decelerates to 2.47 km/s during 30-50 s when the centroidal location passes through the TAPS, but accelerates to 4.64 km/s for 40-70 s, where the centroid is located east of the TAPS (Fig.   6). The latter exceeds the local S-wave velocity of 3.35 km/s at a depth of 5 km where the model plane is set (Table S1). Several studies have proposed that the rupture partially propagated at a supershear speed on the Denali fault, but the estimated velocities varied (Dunham & Archuleta 2004, Ellsworth et al. 2004, Frankel 2004, Asano 2005, Walker & Shearer 2009). Although we cannot directly compare our estimated centroid migration speed with rupture front velocities estimated in other studies , the feature of faster rupture propagation on the eastern ruptured section of the Denali fault is consistent with that resolved by the BP method using teleseismic P-waves (Walker & Shearer 2009).
The focal mechanism projected onto the assumed model plane represents a reverse faulting with a nearly E-W strike direction near the Susitna Glacier fault and a strike-slip faulting with a clockwise rotation of the strike direction from west to east on the Denali-Totschunda fault (Fig. 5d). These results are consistent with surface ruptures observed by field surveys (Eberhart-Phillips et al. 2003) and finitefault models constructed using multiple fault planes (Dreger et al. 2004, Ji et al. 2004, Asano 2005 (Table S2).
In Fig. 8, we show the inversion results with the time-adaptive smoothing constraint. The moment rate function represents a major peak at 26.4 s (Fig. 8a), with the total seismic moment of 1.82 × 10 (& Nm ( 6 8.11). The moment rate function for each basis double-couple component shows that the M3 and M5 dip-slip components and the M1 and M2 strike-slip components are dominant around the peak (Fig. 8b). The time evolution of the potency-rate density along the strike direction (Fig. 8c) shows a large value for 20-40 s, corresponding to the peak of the moment rate function (Fig. 8a). This large slip can be divided into two parts: the vicinity of the rupture front and an area centered about 30 km northeast from the epicenter (Fig. 8c). The rupture front part moves to the northeast as the rupture area expands, but the other part stagnates just before the Xiaoyudong fault (Figs. 8c and d). The snapshot shows that the rupture roughly propagated unilaterally from the epicenter to the northeast, in which reverse faulting was predominant until 60 s from the origin time (Fig. 8d). For 50-60 s, the knot that records the maximum potency-rate density reached about 110 km northeast of the epicenter.
However, for 60-70 s, the location of the maximum knot receded to about 77 km northeast of the epicenter and the moment tensor thereafter turned into strike-slip faulting (Fig. 8d).
The rupture then propagates unilaterally to the northeast direction until the end (Fig. 8d).
The synthesized and observed waveforms are well fitted (Figs. 8e and S6b).
The moment tensor shows predominantly thrust faulting until 60 s from the origin time and then transitions to oblique strike-slip faulting (Fig. 8d). This rake angle rotation toward the northeast from the epicenter is consistent with those estimated by other finite-fault models (Shen et al. 2009, Feng et al. 2010, Nakamura et al. 2010, Fielding et al. 2013, Hartzell et al. 2013. The main rupture, which was propagating northeastward until 20 s, was stagnated in front of the Xiaoyudong fault connecting the Beichuan and Pengguan faults between 20 and 40 s (Fig. 8d), which can be interpreted as the Xiaoyudong fault acting as a geometric barrier preventing the rupture evolution. The timing of the main rupture stagnation is close to the high-frequency emission for 15-25 s near the Xiaoyudong fault detected by the hybrid BP (HBP) method (Okuwaki & Yagi, 2018), which estimates the spatiotemporal distribution of the seismic radiation sources on the assumed fault model by stacking the cross-correlation functions of the observed waveforms and the corresponding Green's functions .
These results are consistent with numerical models and observations that high frequencies are generated when the rupture front velocity changes abruptly (e.g. Madariaga, 1977;Okuwaki et al., 2015).  Fig. 8d). This high-potency-rate density area has significant strike-slip components and thus is distinct from the stagnated rupture with mainly reverse faulting (Fig. S7). On the other hand,  did not identify the significant strike-slip near the rupture front from 20 to 40 s. This suggests that the FFI method is insufficient to estimate a rupture process when the focal mechanism varies. Such a complex rupture process may have been obscured in the conventional FFI methods due to modelling errors caused by the confinement of slip direction and an inappropriate assumption of fault geometry. In contrast, the potency density tensor inversion with high degrees of freedom can successfully estimate the variation of focal mechanisms as well as the complex rupture behavior, which has not been unveiled by the conventional FFI method.

Merits of the time-adaptive smoothing constraint
The numerical experiment highlighted the over-smoothing problem of the unweighted smoothing constraint . Wenchuan earthquake showed that the rupture propagated mainly unilaterally from the epicenter to the northeast, and the obtained focal mechanism changed from thrust to oblique strike-slip faulting, which is consistent with the finite-fault model of previous studies (Shen et al. 2009, Feng et al. 2010, Nakamura et al. 2010, Fielding et al. 2013, Hartzell et al. 2013. Our result clearly distinguished the frontal rupture of strike-slip faulting and the stagnated rupture of reverse faulting near the Xiaoyudong fault from 20 to 40 s, but such a complex rupture process has been difficult to be reproduced with conventional finite-fault models in which the slip vector is fixed to the assumed fault plane. In both applications to real earthquakes, the time-adaptive smoothing constraint successfully estimated complex rupture propagation processes, including fault geometry changes. This is a non-peer reviewed EarthArXiv preprint

23
To further evaluate the effectiveness of the time-adaptive smoothing constraint in the real earthquakes, we performed the potency density tensor inversion with the conventional smoothing constraints  for the 2002 Denali fault earthquake. We then compared these results to those obtained with the time-adaptive smoothing constraint. For these analyses, the parameter settings other than the smoothing constraint were the same as the case with the time-adaptive smoothing constraint. For the moment-tensor-weighted smoothing constraint, " in eqs (6) and (7) was calculated by referring to the GCMT solution (Dziewonski et al. 1981, Ekström et al. 2012 (Figs. 5c and 9b). The presence of this high-potency-rate density area is supported by the location of the sub-event identified from the strong motion data (Eberhart-Phillips et al. 2003, Frankel 2004) and the increase in surface offset from the field survey (Eberhart-Phillips et al. 2003). These results reflect that the unweighted smoothing constraint over-smooths the source signal with particularly large amplitude, as exemplified in the numerical experiment. As a result, the features of the potency-rate density distribution become blurred (Fig. 9b), making the interpretation of the rupture propagation difficult.
The moment-tensor-weighted smoothing constraint weakens the smoothness to the M1 and M2 strike-slip components (Fig. 1) referring to the GCMT solution (Fig. 4).
We show snapshots up to 20 s from the origin time ( Figs. 10a and b), when the dip-slip component is dominant (Fig. 5b) unlike the GCMT solution. The snapshot obtained with the moment-tensor-weighted smoothing constraint shows that high-potency-rate density patches are scattered from 5 to 20 s (Fig. 10a) and the strike direction of the moment tensor rotates nearly 90° from 0-5 s to 5-10 s (inset of Fig. 10a). A high-potency-rate density patch distributed in the southeastern region of the epicenter from 10 to 20 s shows unrealistic focal mechanisms that are sign-flipped from the referenced GCMT solution (Fig. 10a). On the other hand, the snapshot obtained with the time-adaptive smoothing constraint shows a stable potency-rate density distribution while capturing the realistic variation of the focal mechanisms and moment tensor (Fig. 10b). These results show that the moment-tensor-weighted smoothing constraint destabilizes the solution by inappropriate weighting of smoothness constraint when the referenced moment tensor and the focal mechanism for a particular time period deviate significantly, as also shown in the numerical experiment. The moment-tensor-weighted smoothing constraint imposes the same weight over entire source time without distinguishing the slip direction for each basis double-couple component. Therefore, even when the source signal is too small to constrain the focal mechanism, the smoothing can be excessively weak, leading to unrealistic focal mechanisms that are sign-flipped from the referenced moment tensor solution. Such artificial focal mechanism possibly may occur even with the time-adaptive smoothing because the potency density tensor inversion method does not apply the nonnegative slip constraint and the time-adaptive smoothing constraint also weights the five basis double-couple components without distinguishing the slip direction. Nonetheless, since the time-adaptive smoothing constraint controls the hyper-parameters not only in each basis double-couple component but also in the time direction referring to the potency-rate function, it successfully estimates reasonable focal mechanisms even at 10-20 s (Fig. 10b), when the amplitude of the potency-rate function is relatively small (Fig.   5b).
How do we better assign the model planar geometry?
Through the numerical experiment and application to the real earthquakes, we show that the potency density tensor inversion method with the time-adaptive smoothing constraint can be used to project complex rupture propagation onto a horizontal model plane. The potency density tensor inversion method no longer requires an assumption of a rigorous fault geometry, which is a distinct advantage in comparison with the conventional FFI methods; the method provides a robust source model as in the BP method, which also assumes only a horizontal model region to detect a source process. The horizontal model plane, however, cannot properly represent the depth variation of the Green's function because all source knots are located at the same depth . The Green's function of the teleseismic P waveform includes not only the direct phase, but also depth phases reflected at the ground surface near the epicenter, which contain information about the depth of the source. The assumption of a horizontal model plane ignores such source depth information included in the teleseismic waveforms and undermines the high depth resolution, which can be detected by the FFI method using teleseismic waveforms (Yagi 2004, Hartzell et al. 2013). If realistic fault geometry of an earthquake can be approximated by a single model plane, the assumption of an inclined model plane may provide a detailed rupture process including depth resolution. Here, to investigate the effect of the depth variation of the Green's function on the potency density tensor inversion method with the time-adaptive smoothing constraint, we reanalyzed the 2008 Wenchuan earthquake, whose fault geometry is relatively easy to approximate by a single inclined plane, and compared the result with that obtained for a horizontal model plane.
Although finite-fault models with multiple fault planes have been constructed for the 2008 Wenchuan earthquake based on geodetic and geological data (e.g. Shen et al. 2009, Fielding et al. 2013, Hartzell et al. 2013, here we employ a single model plane to simply investigate the effect of depth variation on the Green's function. The assumed model plane has a dip of 32° and an epicentral depth of 16 km, referring to , while fixing the strike to 225° used in the horizontal model plane analysis. The spatial knots of the model plane were reduced from 8.5 to 7.5 km in the depth direction so that the model plane does not extend beyond the ground surface. Because the rupture reached the ground surface (Liu-Zeng et al. 2009, Xu et al. 2009), we did not impose a constraint to make the potency-rate density zero at the top of the model plane. The model settings other than those described above are the same as in the analysis of the horizontal model plane.
In the result of the inclined plane model, the rupture initiating at the hypocenter propagates to the northeast (Fig. 11). The main rupture is arrested in the vicinity of the Xiaoyudong fault zone at about 20 s from the origin time and comes back to the epicenter until about 50 s (Fig. 11). Meanwhile, the near-surface rupture that appears after 20 s continues to propagate northeastward until 60 s as the rupture front expands (Fig. 11). At 60-70 s, a high-potency-rate density area suddenly appears in the center of the Pengguan fault zone and propagates to the northeast until the end (Fig. 11).
The rupture process obtained on the inclined model plane generally agrees with that obtained on the horizontal model plane, suggesting the effect of the depth variation of the Green's function is insignificant, which is consistent with previous studies using the potency density tensor inversion method with the unweighted or the moment-tensorweighted smoothing constraints ). On the other hand, the inclined plane model has the advantage of providing a detailed rupture behavior in the depth direction that is not visible in the horizontal plane model. For example, the rupture extension near the Xiaoyudong fault zone beneath the near-surface rupture around the rupture front at 20-30 s was not clearly captured by the horizontal plane model (Figs. 8d and 11).
We also found minor differences in the rupture behavior not only in the depth direction but also in the strike direction between the horizontal and inclined plane models.
Among them, the most conspicuous difference is the rupture behavior after the stagnation near the Xiaoyudong fault at 20-30 s (Figs. 8d and 11). In the horizontal plane model, the rupture was arrested just before the Xiaoyudong fault until 40 s (Fig. 8d). On the other This is a non-peer reviewed EarthArXiv preprint 27 hand, in the inclined plane model, the rupture that stagnated at about 20 s started to propagate toward the hypocenter until 50 s (Fig. 11). The joint inversion using geodetic, teleseismic and strong motion data estimated bilateral ruptures toward the northeast and southwest directions on the Beichuan fault (Hartzell et al., 2013). The HBP method also showed that high-frequency sources stagnated near the Xiaoyudong fault for 15-25 s and then propagated bilaterally to the northeast and southwest (Okuwaki & Yagi, 2018).
These results support that the back rupture propagation obtained for the inclined plane model is not an artifact. Introduction of the depth variation to the model plane enables us to capture the back rupture propagation.
The success of the inclined plane model shows the importance of including the depth variation to the assumed source knot, which enhances the possibility of capturing a more detailed source process. However, since putting a inclined model plane requires prior information about the fault geometry, it is difficult to apply the model to earthquakes where the fault geometry is poorly constrained due to the lack of surface rupture, such as the 2010 Haiti earthquake (Hayes et al. 2010, Hashimoto et al. 2011).
In addition, a single inclined model plane is difficult to approximate for a fault system with complex geometry, such as the 2002 Denali fault earthquake. Compared to the inclined model plane, the horizontal model plane requires little information about the fault geometry and is widely applicable .
Initial model dependence of the time-adaptive smoothing constraint The time adaptive smoothing constraint is designed so that the relative smoothing becomes weaker in the time period with large amplitude referring to the input potencyrate function, which may lead to an initial model dependence. However, the final solution remained stable even when an unrealistic initial model was used as an input model (Fig.   S4). To examine the effect of the initial selection of an input potency-rate function on the results of the real earthquakes, we tested two types of initial models: uniform slips of the M1 and M5 slip components (Fig. 1), respectively. Comparing the results for both earthquakes, we find no significant difference in the final solution obtained using either of the initial models (Fig. S8). This result shows that the final solution is independent of the initial model parameter settings in the potency density tensor inversion method with the time-adaptive smoothing constraint. Such robustness of the solution suggests that the potency-rate function is stably obtained from the teleseismic body waves that can be approximated by the convolution function of the potency-rate function and the Green's function of centroid.

Conclusion
We developed a potency-density tensor approach of the teleseismic finite-fault inversion,     (Dziewonski et al. 1981, Ekström et al. 2012 and is connected to its centroid position indicated by the cross mark. The rectangle outlines the model plane. The large dots represent surface rupture along the Beichuan, Pengguan and Xiaoyudong faults (Xu et al. 2009), color-coded for each fault. The background topography is derived from the GEBCO 2020 Grid (GEBCO Bathymetric Compilation Group 2020 2020). The inset map shows the regional setting.
This is a non-peer reviewed EarthArXiv preprint