Multiple fault modelling combining seismic and geodetic data: the importance of simultaneous subevent inversions

highlight the strong complementarity of ﬁve different data types: regional seismic data, teleseismic P and S waves, teleseismic surface waves and interferometric synthetic aperture radar (InSAR). The joint data inversions substantially reduce well-known trade-offs between the fault dimensions, fault depth and slip compared to inversions of InSAR or seismic data alone. The synthetic inversions show that local/regional seismic data are key to constrain the lowest magnitude subevents, notably for the two-fault conﬁguration with downdip listricity, which cannot be well constrained by the InSAR or teleseismic data alone. This suggests that downdip segmentation can be difﬁcult to detect in the absence of local/regional seismic data.


I N T RO D U C T I O N
While some moderate to major earthquakes appear to be relatively simple, occurring on a single fault plane, there is increasing evidence of complex earthquake behaviour. For instance, the M w 7.3 1992, Landers earthquake (e.g. Fialko 2004), the M w 7.1 1999, Hector Mine earthquake (e.g., Hauksson et al. 2002) or the 2010-2011 seismic sequence in Canterbury, New Zealand (e.g. Atzori et al. 2012;Elliott et al. 2012) are just some examples of events that involved the rupture of multiple faults. Some earthquakes come as a surprise, rupturing previously unmapped faults (e.g. Berberian 1979;Geersen et al. 2015), and surface geology observations can give a limited view of the source process (Talebian et al. 2004;Fialko et al. 2005). Geophysical techniques such as seismology and geodesy are thus needed for detailed investigations of the rupture, the seismic cycle and seismic hazard.
One of the simplest descriptions of an earthquake is a point source characterized by a space-time centroid and a moment tensor. Point source models are routinely estimated for earthquakes on a global scale (e.g. Dziewoński et al. 1981;Ekström et al. 2012;Duputel et al. 2012b). However, more advanced source descriptions are needed to understand the detailed source process and the physics of earthquakes.

Seismo-geodetic multiple fault modelling
Since the 1980s, several methods have been developed to constrain the space-time distribution of slip along a fault plane (e.g. Kikuchi & Kanamori 1982;Das & Kostrov 1990;Ji et al. 2001Ji et al. , 2002. Under a number of assumptions (e.g. constant rise-time and rupture speed, simple and identical elementary slip function for all points on the fault), one can use the representation theorem to perform linearized inversions for slip, such as in classical multipletime-window approaches (e.g. Olson & Apsel 1982;Hartzell & Heaton 1983). However, the assumptions used in these approaches can be too strong and uncertainty quantification is typically limited because the parameter space is not extensively explored. On the other hand, fully non-linear distributed slip inversions (e.g. Ji et al. 2002;Liu & Archuleta 2004) have the advantage of simultaneously inverting for all the distributed slip source parameters including, for example the rise time and rupture speed. Recent developments include the use of Bayesian approaches fully exploring the parameter space (e.g. Monelli & Mai 2008;Monelli et al. 2009;Razafindrakoto & Mai 2014;Minson et al. 2014a;Duputel et al. 2015;Duputel & Rivera 2017). Moreover, the United States Geological Survey (USGS) now routinely publishes finite fault distributed slip maps for large magnitude events based on the technique of Ji et al. (2002). Despite these efforts, a key weakness of distributed slip inversions is their non-uniqueness, particularly when only teleseismic data are used in the inversions (Das & Kostrov 1990;López-Comino et al. 2015). For example, distributed slip maps produced by different research groups for a given event can be quite different, even when the same data sets are used (e.g. Weston et al. 2012;Mai & Thingbaijam 2014;Mai et al. 2016). Near-field data from strong-motion sensors can be beneficial for modelling the spacetime evolution of seismic slip on fault planes (e.g. Olson & Apsel 1982;Hartzell & Heaton 1983), and can contribute to reducing the non-uniqueness of finite-fault distributed slip maps in inversions based only on far-field data. However, the availability of near-field stations is often limited.
In order to overcome the non-uniqueness and stability issues associated with distributed slip models, so-called intermediate techniques have been developed. For example, Vallée & Bouchon (2004) developed a method using elliptical slip patches to model earthquake sources, which was recently revisited by other research groups (Di Carli et al. 2010;Ruiz & Madariaga 2011;Ulrich & Aochi 2015). McGuire et al. (2001) used second-degree polynomial moments of earthquake space-time distributions to estimate earthquake fault dimensions, location, duration as well as the rupture velocity. Another way to model an earthquake is by using multiple point sources rather than a single point source (e.g. Kikuchi & Kanamori 1991;Tsai et al. 2005;Duputel et al. 2012a;Duputel & Rivera 2017). Such techniques have received increasing attention in the past decade, but they tend to be limited to long-period seismic data and thus have been restricted mostly to the study of large earthquakes. Additional data types, including shorter period data, are needed to investigate multiple fault rupture processes for moderate magnitude events. Moreover, the effect of 3-D earth structure on multiple point source inversions is starting to receive some attention (e.g. Duputel & Rivera 2017).
The M w 7.2, 1992 Landers earthquake marked the beginning of the use of space-born geodetic data such as GNSS and InSAR in earthquake source modelling (e.g. Massonnet et al. 1993). In the past decades there has been an explosion in the use of various data types such as regional seismic waveforms, strong motion data, teleseismic waveforms, GPS, InSAR and other geodetic data (Konca et al. 2010;Salichon et al. 2003;Delouis et al. 2010;Funning et al. 2014). Combining different types of data is important to reduce the non-uniqueness of the source inversions and helps to better constrain the earthquake source models (Weston et al. 2014). For example, while seismic data have high temporal resolution, InSAR and GPS data have high spatial resolution, highlighting the complementarity of these data types (e.g. Weston et al. 2012).
Different sources of data can be used to model earthquake source parameters in a multistage inversion whereby, for example, a fault with static slip is first modelled only with geodetic data and in the following step the temporal slip history is modelled with seismic data. Such multistage inversions are common in the literature (e.g. Hernandez et al. 1999;Semmane et al. 2005;Custódio et al. 2009), there are disadvantages in these approaches such as not accounting for the uncertainty of the parameters from the inversion of single data sets. Possible trade-offs between fault parameters and data sets may be hidden with such approaches. On the other hand, an advantage of multistage inversions would be the reduced computational cost compared to a simultaneous inversion of different data sets.
In geodesy, in addition to distributed slip inversions where the fault geometry is often fixed, several studies have used InSAR and/or GPS data to determine multisegment complex fault geometries where each segment may have different orientations and average slips (e.g. Wright et al. 1999Wright et al. , 2001Xu et al. 2010;Atzori et al. 2012;Kobayashi et al. 2013;Akoglu et al. 2018). Such studies often use geological information to reduce the number of parameters determined (e.g. Feng et al. 2010;Fielding et al. 2013;Hamling et al. 2017), and iterative inversions are typically used (e.g. Wright et al. 1999;Nishimura et al. 2008;Huang et al. 2016). Specifically, first the inversion for single fault parameters (e.g. strike, dip, rake, spatial location, fault dimensions and average slip) is carried out. If the misfit between the observations and the simulations is large, an additional fault is added while keeping the parameters of the first fault fixed to the solutions found in the first iteration. More faults are then iteratively added until the data misfit does not decrease significantly.
The multiple fault modelling technique presented in this study goes beyond previous inversion schemes determining average multifault parameters by incorporating a wider range of data types than before (local, regional and teleseismic data, as well as InSAR and GPS) and by incorporating a 3-D earth structure model in the forward modelling of the seismic data. The technique enables the inversion for parameters of all faults at the same time (which we refer as 'simultaneous inversion' throughout this manuscript). Moreover, it is also well suited for inversions where one fault is modelled after the other (which we refer as 'iterative inversion', as described in the previous paragraph).
Our approach is data-driven, whereby we successively add complexity to the source models as required by the data. We test for the first time whether the widely used iterative inversion approach is unbiased with respect to the final model. Additionally, we evaluate the effect of data noise on the models and the associated uncertainties, and the different constraints provided by the different input data types.

Model formulation
The new multiple fault modelling technique builds on the joint inversion scheme introduced by Funning (2005) and further developed by Weston et al. (2014) to create single fault earthquake 960 M. Frietsch et al. source representations using seismic and geodetic data. While Weston et al. (2014) inverted for a single double-couple source with nine source parameters-fault strike, dip, rake, average slip, length, width, centroid longitude, latitude and depth-here two additional source parameters are determined: the time-shift to the centroid time and the compensated-linear-vector-dipole (CLVD) component, following the approach of Tape & Tape (2015). Furthermore, the new scheme includes the ability to invert for the parameters of multiple faults simultaneously, as well as the traditional iterative approach often used in source inversions, as explained in the introduction (e.g. Wright et al. 1999;Nishimura et al. 2008;Huang et al. 2016). The algorithm can also model not only a single earthquake via multiple faults, but also multiple earthquakes captured by a single InSAR image (e.g. foreshocks and the main shock), with the seismic data being used to discriminate between the various events.
The search for the source parameters is bound-constrained and uses the Powell algorithm (Powell 1964) as a local optimization scheme with multiple Monte Carlo restarts in order to find the global minimum of the L 2 -norm misfit (see Fig. 1 with a flowchart summarizing the source inversion algorithm). The misfit function m 2 is based on the difference between the theoretical seismograms t S and the observed seismograms d S , as well as the difference between the theoretical geodetic t G and observed geodetic d G displacements (with the weighting factors for the seismic α S and geodetic α G components): Considering a seismic moment tensor f : the theoretical seismograms t S can be written as: where K is a matrix with the six sensitivity kernels of the seismic waveforms with respect to each component of the moment tensor for a given space-time centroid (e.g. Ferreira & Woodhouse 2006). The kernels are the partial derivatives of the synthetic waveforms with respect to the moment tensor components: The seismic forward modelling scheme is flexible; the theoretical seismograms t S and the seismic kernels K can be computed with any suitable waveform modelling technique. Here, we use the spectral element method to model seismic waves propagation in a 3-D earth model with the package SPECFEM3D GLOBE (Komatitsch & Tromp 1999, 2002a, which is a highly accurate technique for seismic forward modelling in realistic 3-D earth models. Given the linear relationship between the moment tensor and the seismic waveforms (eq. 3), the sensitivity kernels can be pre-computed and stored before being used in the source inversions. The geodetic forward modelling uses the elastic dislocation theory for the displacement of a rectangular fault in a homogeneous half-space (Okada 1985). Hence, there is a difference between the complexity of the earth structure model used in the seismic and in the geodetic modelling (3-D model versus half-space). This issue will be discussed in Section 7.

Inversion scheme
As a first step, we invert the different seismic and geodetic data sets separately to explore their sensitivity to the different source parameters. Then, the weighting factors (α S and α G ) in the inversions with multiple data sets are adjusted so that each data set contributes equally to the misfit function m 2 , by taking the inverse of the misfit values from the separate inversions as a weighting factor. The weighting factors need to be adjusted for every earthquake. For extreme cases such as very noisy data or a scarce data set, this needs to be taken into account for the determination of sensible weighting factors and in certain circumstances data sets might even need to be removed. A more rigorous weighting needs to be applied in future work, as discussed in Section 7. The simultaneous, multiple fault inversion algorithm successively adds subfaults to the modelling following a data-driven approach. Every source inversion starts by considering a single fault. The number of subfaults is then increased one by one until the improvement in data misfit is smaller than 5 per cent or if the seismic moment of the additional subevent becomes smaller than the previous one by an order of magnitude or more. The 5 per cent threshold in the data misfit criterion is consistent with previous seismic studies that used similar thresholds taking data noise into consideration (e.g. Debayle & Ricard 2012;Parisi et al. 2014).
In principle, the algorithm works for an arbitrary number of subfaults, and several configurations of two-and three-fault models were tested. However, for conciseness, the remainder of this study is focused on two-fault examples. The number of restarts required by the algorithm depends on the complexity of the source model. Having extensively experimented with the number of Monte Carlo starting points, we found that 200 to 1000 Monte Carlo restarts are enough for a single fault model to converge, while a twofault model needs 100 000 different Monte Carlo starting points. We found that in order to ensure that convergence is achieved, the number of required Monte Carlo restarts depends strongly on the number of fault parameters. Our tests showed that when doubling the number of parameters, the number of restarts needs to be squared.

S Y N T H E T I C C A S E S T U D I E S
In this study we perform synthetic inversion tests for two input artificial earthquakes, each comprising two subfaults with: (i) downdip listricity, a geometry inspired by the spatial distribution of the aftershocks of the M w 6.0, 21 February 2008 Wells, Nevada event in the Basin and Range province of the western USA (Smith et al. 2011) and (ii) along-strike branching broadly compatible with surface observations of the M w 7.1, 16 October 1999 Hector Mine earthquake in California (Simons et al. 2002). Fig. 2 shows the earthquake scenarios considered. Since they are based on information from real events, the artificial ruptures considered in this study exhibit realistic source complexity. The seismic and InSAR data configurations considered for the two events are based on those for the M w 6.0, 2008 Wells event, which was well captured by high-quality In-SAR observations as well as by extensive local and teleseismic data (Nealy et al. 2017). We compute synthetic data for these scenarios as explained in Section 5.

DATA
The inversion scheme described above can incorporate various different data types, such as geodetic data (e.g. GNSS displacements,  InSAR, levelling, Lidar), geological field data (e.g. surface fault offsets and/or mapped traces) and seismic data (e.g. strong motion, GNSS seismology, local, regional and teleseismic waveform data). This study focuses on noise-perturbed synthetic tests with InSAR, regional seismic and teleseismic data sets computed for the two earthquake scenarios considered. We note that in this study we refer to regional seismic data as waveforms recorded at epicentral distances between 0 • and 10 • . Moreover, we consider teleseismic data in the epicentral distance range of 30 • and 140 • .

InSAR data
We generate synthetic interferograms with acquisition geometries based upon the Envisat satellite in both ascending and descending orbits. We use the analytical solution of Okada (1985) to estimate the 3-D displacement field for each synthetic event and project them into the satellite's line of sight. We add realistic noise to the synthetic InSAR data by convolving a random vector of specified amplitude with the 1-D covariance relationship estimated from an undeformed portion of the real interferograms of the M w 6.0, 2008 Wells earthquake (e.g. Hanssen 2001;Wright et al. 2003). Fig. 3 shows the noise-free synthetic interferograms used compared to noise-perturbed examples, as well as the L 2 -norm misfit between the two.

Seismic data
The distribution of stations for all seismic data sets corresponds to real stations from the global seismic network and from the temporary USArray, which was deployed in the region of the M w 6.0, 2008 Wells earthquake at the time of the event. SPECFEM3D GLOBE (Komatitsch & Tromp 1999, 2002a is used to compute synthetic waveforms, taking SGLOBE-rani (Chang et al. 2015) as the mantle model and CRUST2.0 (Bassin et al. 2000) as the global crustal model. Those 3-D earth structure models for seismic data are in contrast to the homogeneous half-space used for the geodetic forward modelling. This difference in earth structure should not have a severe impact on the overall conclusions of this study and is discussed in Section 7. We perform two different seismic simulations in order to save computing time: a global simulation of 2-hr-longteleseismic waveforms (with epicentral distance between 30 • and 140 • ) accurate down to a wave period of T ≈ 17 s, and a regional simulation of 30-min-long waveforms (with epicentral distance between 0 • and 10 • ) accurate down to T ≈ 10 s.
We compute three-component teleseismic waveforms for the locations of stations from the GSN (IU, II), Caribbean (CU), China (IC), Geoscope (G), Geofon (GE) and MedNet (MN) networks. Likewise, regional synthetic seismic waveforms are computed for Transportable Array stations of the USArray in the Western USA, in the range between 31.3-49.1 • N and 106.4-123.6 • W. The distribution of stations for regional seismic data, teleseismic P waves, 962 M. Frietsch et al. Figure 2. Two-fault geometry and source parameters used in the synthetic tests performed in this study for: (a) downdip listricity, whereby the largest magnitude subevent (blue) has M w 6.0 and the lowest magnitude subevent (red) has M w 5.3 and (b) along-strike branching, whereby the largest magnitude subevent (blue) has M w 6.8 and the lowest magnitude subevent (red) has M w 6.6. The two subfaults are illustrated by dark red and blue areas, the lines in corresponding colours mark the projected surface breakout and the arrows show the slip direction. The beach-balls correspond to each subevent according to their colour. The source inversions determine eleven parameters for each fault: spatial centroid location (latitude, longitude, depth; yellow circles), fault width, length, strike, dip, rake, average slip, time-shift and the CLVD component. teleseismic S waves, and teleseismic surface waves is shown in Fig. 4. Table 1 shows the seismic processing parameters used, such as filter range, window length, epicentral distance range and number of seismic traces.
We create realistically noisy seismic data by adding real waveform noise to our synthetic seismograms. We cut a window from background noise recorded before the P-wave arrival in seismograms of the M w 6.0, 2008 Wells earthquake. In order to ensure that the levels of noise used are realistic, the noise is normalized and scaled by a specified amplitude to achieve a broad range of signal-to-noise ratios (see Fig. S2b), and is consequently added to the synthetic waveforms. Noise-free versus noisy synthetic waveforms for the along-strike branching event based on the M w 7.1, 1999 Hector Mine earthquake are compared in Fig. 5 (see Figs S4-S6 in the supplementary materials for waveform comparisons between noisy and noise-free synthetics for the source model with downdip listricity). Using different earth structure models in the forward modelling schemes for geodetic (homogeneous half-space) and seismic data (3-D model) should not affect the main conclusions of this study and is discussed in Section 7.

R E S U LT S
This section presents the results from synthetic inversion tests designed to evaluate the performance of the new source inversion algorithm, based on the fault configurations shown in Fig. 2. We show first a set of tests focused on the robustness of the source inversions in the presence of data noise. Secondly, we investigate the limits of the method for detecting lower magnitude subevents, by varying the time-shift between the two events. Finally, we evaluate the differences in source parameters obtained from iterative versus simultaneous source inversions.

Effect of data noise
Here we investigate the capability of the simultaneous inversion scheme to recover the input fault model in the presence of noise. In particular, we examine the resolving power of each input data set. We invert individually and jointly 100 noise-perturbed data sets for each synthetic case study. The amplitude of the noise in each realization is drawn from a Gaussian distribution with a standard deviation in the range of 3.4-5.9. This leads to a nearly Gaussian distribution of misfits between the 100 noise-perturbed and noise-free synthetic InSAR displacements for both ascending and descending tracks (Fig. S2a). Likewise, the 100 sets of noise-perturbed synthetic waveforms have nearly Gaussian signal-to-noise ratio distributions with a median of 7 (Fig. S2b). The overall level of noise used in the synthetic tests is chosen such that the signal-to-noise ratios achieved are comparable to those in real InSAR and seismic data.
We invert simultaneously for the 22 source parameters of the two earthquake scenarios considered (Fig. 2) using the inversion scheme presented in Section 3. We carry out both separate and joint inversions of the various 100 noise-perturbed synthetic data sets, that is, combined seismic (regional seismic, teleseismic P waves, teleseismic S waves, teleseismic surface waves), and combined InSAR (ascending and descending interferograms). Fig. S10 shows the source parameters obtained from the inversion of the combined seismic data set for each subfault (hereafter the 'seismic-only' inversion) for the two-fault model with downdip listricity. The results for the 100 noise-perturbed data sets are shown as trade-off plots, that is, scatterplots between each pair of source parameters, and histograms of the source parameter distributions (e.g. Wright et al. 2003;Weston et al. 2014). We identify trade-offs between source parameters as positive or negative correlations in individual scatterplots; for example, Fig. S10(bottom panel) shows a clear correlation between width and slip.
The larger (M w 6.0) of the two subevents is well constrained (Fig. S10, top). Parameter trade-off plots and histograms for each individual fault parameter show tight clusters of solutions around the input value in each case, suggesting that the method is capable of retrieving these well, with the exception of the CLVD component, which has a wide range of retrieved values and shows strong trade-offs with the other source parameters. The CLVD component seems to be an 'error absorbing' parameter that generally increases  . Station distribution for the seismic data sets used in this study: (a) regional seismic waves; (b) teleseismic P waves; (c) teleseismic S waves; and, (d) teleseismic surface waves. The diagrams show the source-receiver azimuth range from 0 • to 360 • ; while the radial direction shows the epicentral distance range in degrees. with higher levels of noise, while the other parameters are relatively stable. The smaller subevent is also well constrained with trade-offs only evident for strike, dip and CVLD component, as well as some uncertainties (of about 15 • for fault strike, 10 • for dip, 10 • for rake and 98 per cent ( ν ≈ 0.33) for the CLVD ν component; Fig. S10, bottom). As expected, the seismic data have high temporal resolution, notably the higher frequency regional seismic waveforms, and are very sensitive to the time-shift between the two subfaults. This is evident from inversions carried out using each data set separately (see Section 3), which are always performed before the joint inversion of all data sets. We present the results of the 'InSAR-only' inversions in Fig. S11. The distributions of source parameters are broader than in the seismic-only inversion and reveal larger fault parameter uncertainties for the shallower, smaller subevent (Fig. S11, bottom panel) than for the deeper, larger subevent (Fig. S11, top panel). There are substantial trade-offs between all source parameters within each subevent, and also clear trade-offs between the parameters of the two subevents. The InSAR data alone cannot fully retrieve the source parameters for this fault configuration in the presence of noise. The problem seems to be strongly non-unique, as the misfit varies only slightly for a large range of different earthquake models. The non-uniqueness may be due to the decrease in resolution with depth for InSAR data as well as the correlated noise which could change the shape of the main deformation lobe and therefore significantly affect the estimated dip angle for dip-slip faults. The joint seismic/InSAR inversion leads to a significant reduction of uncertainty in the fault parameters (Fig. 6) compared to the individual data inversions. The trade-off plots for the joint model show much more clustered results, with sharp peaks in the histograms, especially for the larger, deeper subevent (Fig. 6, top  panel). The input fault parameters are overall well retrieved; minor trade-offs remain in the fault slip, width, depth and length, as seen in previous studies (Weston et al. 2014), as well as in the CLVD component. These parameters exhibit larger errors than the others (≈0.3 m in fault slip, ≈9 km in fault width, ≈4 km in depth and ≈8 km in fault length), but they nevertheless are smaller than those obtained in the single data inversions (Figs S10 and S11). The shallower, smaller subevent (Fig. 6,  the joint data inversions compared to the InSAR-only solutions. Despite having a high artificial CLVD component and suffering from the same parameter trade-offs as the larger magnitude subevent, the other source parameters are well retrieved. When considering the two-fault source model with along-strike branching, we find that overall the various source parameters for the two subevents are well recovered in the joint data inversions (Figs S14) and in the InSAR data inversions (Fig. S13). The seismic-only inversions also retrieve the input parameters well (Fig. S12).

Comparison of iterative and simultaneous inversion of multiple faults
In this section we perform synthetic inversion tests for the two-fault earthquake scenarios shown in Fig. 2 using: (i) an iterative approach as typically used in previous studies (see Section 1 for details), whereby only 11 source parameters are determined in each iteration and (ii) the simultaneous inversion for two faults implemented in this study, where 22 source parameters are jointly determined.
While in the previous section we used 100 noise-perturbed data sets, in this section we use a single synthetic data set for each data type (i.e. combined seismic and combined InSAR) whose noise amplitudes correspond to the median of the distributions shown in Fig. S2. Table 2 compares the results obtained from an InSARonly inversion with those from the joint inversion of all data types, for the source model with downdip listricity (Table S1 shows the corresponding results for the along-strike branching case study). The iterative InSAR inversion overall retrieves well the input fault parameters for the larger magnitude subevent (Table 2). However, there are substantial differences for the lower magnitude subevent compared to the input values, such as a difference of ≈35 • in fault strike, ≈54 • in dip, ≈12 • in rake, 0.1 in moment magnitude and ≈15 km in centroid location. A similar behaviour can be observed in the results from the iterative joint data inversion. While the parameters of the larger magnitude subevent are well retrieved, there are substantial discrepancies for the lower magnitude subevent (≈22 • for fault strike, ≈2 • for fault dip, ≈41 • for fault rake and M w = 0.1 for moment magnitude. As expected, the differences from input values in the joint inversion are smaller than in the InSARonly inversion, but overall the iterative inversion approach leads to an incorrect solution for the lower magnitude subevent. Similar parameter trade-offs are found to those discussed in the previous section, notably between between fault slip, width, length and depth. In order to examine the variability of solutions obtained from the 966 M. Frietsch et al.  Figure 6. Trade-off plot for 100 synthetic test inversions for the source model with downdip listricity using joint seismic and InSAR data with different levels of noise added. The circles mark the respective best-fitting model. The yellow star marks the inversion's solution with the lowest misfit. The light blue triangles are the input parameters for the subevent with the larger magnitude whereas the pink triangles are the input parameters for the subevent with the lower magnitude (see Table 2); the same colour code applies for the vertical dashed lines in the histograms. We plot the following parameters: the fault strike, dip, rake, slip, the centroid longitude and latitude (for UTM zone 11T), the fault length, width, the centroid depth of the fault, the time-shift and the CLVD component. The beach-ball of the input model is shown for each subevent. Table 2. Comparison of the fault parameters for iterative and simultaneous two-fault inversions for the input two-fault model with downdip listricity. This is a synthetic test corresponding to the two-fault earthquake scenario shown in Fig. 2(a), with noise added for the four seismic data sets (teleseismic P waves, teleseismic S waves, teleseismic surface waves and regional waveforms) and the InSAR data sets. The fault parameters are the centroid longitude and latitude (for UTM zone 11T), the centroid depth, the fault length, width, slip, the fault strike, dip, rake, the double couple percentage, the time-shift to the GCMT time, the L 2 -norm misfit, the seismic moment M 0 (Nm), the magnitude and the corresponding beach-ball (Stereo). The colours of the beach-balls correspond to the larger magnitude subevent (blue) and the lower magnitude subevent (red iterative inversions that fit the data well, Fig. 7 shows the best-fitting solutions not exceeding 20 per cent of the best misfit. For the larger magnitude subevent (Fig. 7, top panel), the lowest misfit solutions tend to cluster around the correct input values, apart from the solutions for slip, which show more variability. On the other hand, for the lower magnitude subevent (Fig. 7, bottom), the solutions are distributed on the whole parameter space for the majority of the fault parameters, displaying substantial errors and uncertainties. The results from the simultaneous inversions based on InSAR data alone (Table 2) underline the difficulties in using noisy data to model the source model with downdip listricity. Similar to the iterative inversions, the retrieved fault parameters differ significantly from the input model for both subfaults and the two-fault geometry considered is not properly recovered. The solution for the larger magnitude subevent shows differences of 2.3 km in fault length, 1.9 km in fault width and 8.6 • in fault strike, while the other parameters are retrieved reasonably well. The solution obtained for the lower magnitude subevent overestimates the moment magnitude by M w = 0.2, which results from the systematic error in fault width of 1.3 km and in fault slip of 0.04 m (27 per cent higher). The fault dip is 22.1 • shallower than in the input model and the fault rake has an systematic error of 18.8 • . On the other hand, as shown in the previous section, the simultaneous joint source inversion recovers well the two subfault solutions. All the parameters are well retrieved, with only some slight systematic errors (e.g. an error of ≈12 • in fault rake for the lower magnitude subevent) that are within previously reported expected error bounds in source parameters (e.g. Weston et al. 2012). Table S1 shows that the iterative inversions for the source model with along-strike branching lead to smaller errors than for the event with downdip listricity. Most source parameters are well recovered, with the largest systematic errors being of 0.1 for moment magnitude, 22 per cent in non-double couple component for the smaller magnitude event and ≈1 m in slip.

Effect of the interevent time between the (sub)events
In order to determine the limits of detectability of the lower magnitude subevent considered in the source model with the alongstrike branching event (see the corresponding fault configuration in Fig. 2b), we carried out source inversions for a series of different input interevent times between 0 and 100 s. Since the regional seismic data are key to constrain the subevent, we perform inversions only using regional seismic data with a signal-to-noise ratio of 7. Fig. 8 shows the differences between input and output parameters obtained as a function of the interevent time in the iterative inversions. It shows that for interevent times larger than about the dominant wave period of the regional seismic data (T ≈ 25 s) all the parameters are well constrained. Smaller subevent intervals lead to the following average sytematic erros in the source parameters: ≈10-25 • in fault strike, ≈5-10 • in fault dip, ≈10-15 • in fault strike, ≈5-20 per cent in CLVD component, ≈5 × 10 18 , and ≈ 3 s in the interevent time (with an outlier of 21 s time difference possibly due to a cycle skip). The time difference from which the iterative inversion technique works well is larger than the expected duration of an earthquake with the magnitude considered. Thus, the iterative approach is valid when there is a large temporal separation between the two subevents. In contrast, the simultaneous inversions (Fig. S1) recover well the input parameters for much smaller temporal separations in order of the sampling rate.

D I S C U S S I O N
This study shows that joint inversions of seismic and geodetic data in the presence of noise constrain two-fault models reasonably well when inverting for all source parameters simultaneously. Spatial source parameters such as centroid latitude and longitude are well retrieved in InSAR and joint data inversions (Fig. S11), being a strength of InSAR with its high spatial resolution. The fault geometry and rupture orientation are also fairly well constrained by InSAR, seismic and joint data inversions. The time-shift is constrained by the seismic data (especially the regional data set) and is well determined in both seismic and joint data inversions. We find that the CLVD component is difficult to constrain in all cases and shows errors that can go up to 98 per cent ( ν ≈ 0.33). Our work suggests that the CLVD component mostly acts as an 'error absorbing' parameter in the inversions, which is very sensitive to data noise. The CLVD component quantifies the deviation from a pure double-couple source, and indicates either source complexity, such as an earthquake rupturing on multiple faults with different fault geometries, non-planar faults, fluids or anisotropy in the source volume Miller et al. 1998). Alternatively, the CLVD component may be spurious due to inaccurate earth structure models, or imprecise waveform modelling techniques (Henry et al. 2002;Ferreira et al. 2011). The large CLVD components retrieved in the synthetic inversion tests performed in this study are due to the effect of noisy data on the inversions. This suggests that when performing multiple fault inversions of real earthquakes it may be difficult to distinguish whether a large CLVD component is real or spurious. The InSAR modelling, in contrast to the seismic part, does not include the CLVD component, which should be the subject of future work. We do not consider a tensile mode rupture in this study, although the displacement modelling method of Okada (1985) would allow to account for a fault opening with slip perpendicular to the fault plane, in addition to a shear failure with slip in the fault plane. This shear-tensile-crack (or sometimes called 'crack + double-couple') source description (e.g. Dufumier & Rivera 1997;Minson et al. 2007;Tape & Tape 2012 would affect both the isotropic and CLVD component of the seismic moment tensor (e.g. Vavryčuk 2001Vavryčuk , 2011. Therefore, changing the seismic parametrization from a deviatoric source to a shear-tensile-crack type source as in Minson et al. (2007), where seismic and geodetic data is modelled with a more consistent source description, would be beneficial in future work. Another feasible way would be to change the source parametrization in the seismic fault modelling to a full moment tensor with a decomposition to a shear-tensile-crack source (e.g. Tape & Tape 2013). Likewise, with the change of the geodetic forward modelling scheme to, for example Zhu & Rivera (2002), both seismic and geodetic source parametrizations could be based on the same full moment tensor where a deviatoric constraint can be easily included in the inversion scheme. All three options would open the possibility to study events with a tensile mode (e.g. explosions, collapses, volcanic events, as well as induced and natural geothermal events) in a more sophisticated way. As expected, when considering the source model with downdip listricity (Fig. 2a), the subevent with the lower magnitude (one magnitude unit lower than the other subevent) is harder to constrain. In particular, when using InSAR data alone in the inversions, the corresponding source parameters cannot be reliably retrieved. In addition, substantial tradeoffs between the two subfault parameters are observed. This study finds that regional seismic data are important for capturing and constraining lower magnitude subevents in faulting scenarios with downdip segmentation. Thus, this study suggests that in the absence Seismo-geodetic multiple fault modelling . Trade-off plots for the best-fitting solutions not exceeding 20 per cent of the best misfit from the synthetic test inversions for the source model with downdip listricity using joint seismic and InSAR data with median levels of noise added (see main text for details). The grey circles represent the solutions. The yellow star marks the inversion's solution with the lowest misfit. The light blue triangles are the input parameters for the subevent with the larger magnitude whereas the pink triangles are the input parameters for the subevent with the lower magnitude (see Table 2); the same colour code applies for the vertical dashed lines in the histograms. We plot the following parameters: the fault strike, dip, rake, slip, the centroid longitude and latitude (for UTM zone 11T), the fault length, width, the centroid depth of the fault, the time-shift and the CLVD component. The beach-ball of the input model is shown for each subevent. Figure 8. Differences between the input and output fault parameters from synthetic iterative inversion tests whereby the time difference between the two subevents is successively varied from 0 to 100 s for the along-strike branching event (see the corresponding fault configuration in Fig. 2b). In this synthetic test, the two fault sources have an input time-shifted between 0 and 100 s. The differences are shown for fault strike (a), fault dip (b), fault rake (c), the CLVD component (d), seismic moment (e) and the interevent time observed (f). The colours correspond to the larger magnitude subevent (blue) and the lower magnitude subevent (red). of such local/regional data, downdip segmentation may be difficult to detect. Yet, such segmentation has been suggested in numerous regions worldwide (e.g. Iran, California, Italy, China, Turkey; Elliott et al. 2011;Bent & Helmberger 1989;Ross et al. 2017;Chiaraluce et al. 2011;Wang et al. 2009;Cakir & Akoglu 2008); it can lead to very distinct seismic hazard scenarios (Ofoegbu & Ferrill 1998;Boncio et al. 2004;Passone & Mai 2017) such as the decrease of peak ground velocities on the footwall and the increase on the hanging wall with increasing listricity. Constructive interference of seismic waves due to a listric fault can lead to over twice higher peak ground velocities than those observed for a planar fault. Downdip segmentation has important implications in terms of tectonics and earthquake physics (Ofoegbu & Ferrill 1998;Oglesby & Archuleta 2003). It can change the stress field and could therefore influence rupture complexity and seismic slip in future earthquakes. Our inversions show that the source model with along-strike branching (Fig. 2b) is better constrained than the model with downdip listricity. We show that harnessing the spatial resolution of InSAR in the joint data inversions is key to enhance the robustness and accuracy of the estimated source parameters of these events compared to single data inversions. Hence, this work shows that in addition to the relatively simple single fault models considered by Weston et al. (2014), the study of along-strike branching events greatly benefits from using multiple seismic and geodetic data sets with complementary sensitivity. In all InSAR and joint data inversions, we observe a trade-off between fault slip, width and the depth of the fault; this is a well-known issue in inversions with InSAR data (e.g., Simons et al. 2002;Sudhaus & Jónsson 2009;Weston et al. 2014). Nevertheless, compared to single data inversions, the joint inversions lead to narrower peaks (i.e. smaller uncertainties) in the parameter distributions and thus to more robust solutions. Moreover, combining the various data sets in the inversions helps reduce parameter tradeoffs, such as between fault dip, slip and length. Whenever reliable a priori information is available (e.g. surface geology constraints on fault strike and slip), it may also help further stabilize the inversions and to reduce parameter trade-offs (Walters et al. 2018). Since many studies use iterative source inversion approaches both using geodetic (e.g., Wright et al. 2001;Nishimura et al. 2008;Huang et al. 2016) and seismic data (e.g. Kikuchi & Kanamori 1991;Sokos & Zahradník 2008;Zahradník & Sokos 2013;Quintero et al. 2014), we compare the results from iterative and simultaneous multiple fault inversions. This study shows the need to invert simultaneously for all subfault parameters of the source model with downdip listricity whereby one of the subevents has a moment magnitude ≈0.7 lower than the other, since the iterative approach leads to significant errors in the estimated source parameters (e.g. for InSAR-only modelling, errors of ≈35 • in fault strike, ≈54 • in dip, ≈12 • in rake and 0.1 in moment magnitude; and for the joint data modelling, errors of ≈22 • for fault strike, ≈2 • for dip, ≈41 • for rake and 0.1 for moment magnitude).
A limitation of our geodetic modelling scheme is the assumption of an isotropic, homogeneous, elastic half-space (Okada 1985). Previous comparisons between the effects of layered versus homogeneous half-spaces on the geodetic modelling showed that the latter may lead to systematically shallower depths, but with a bias not exceeding 30 per cent (Savage 1987;Lohman 2002;Fialko 2004;Hearn & Bürgmann 2005), which depends on the source mechanism, and may be smaller than the effects, for example, of data errors. We expect that a depth biased by 30 per cent has a negligible effect on the seismic data as even for regional waveforms this systematic error represents only a small fraction of the relatively long wavelengths used in this study. Nevertheless, for internal consistency, in the future we will use a layered elastic structure in our modelling of the geodetic data using the same local layered structure in the source region as in the 3-D earth models used in the seismic modelling (e.g. Zhu & Rivera 2002). In addition, a fully local 3-D earth model may be required for accurate modelling of the Green's functions in areas with substantial lateral variations in elastic structure (e.g. Barbot et al. 2009).
The rectangular shape of the dislocation source used (Okada 1985) is in general different from observed fault structures in the field (e.g. Fossen & Rotevatn 2016;Mildon et al. 2016;Civico et al. 2018). Geodetic deformation patterns for shallow faults (e.g. rupturing the surface) show evidence for non-rectangular dislocation sources (e.g. Maerten et al. 2005;Furuya & Yasuda 2011). Walters et al. (e.g. 2009) see an improved prediction of the surface rupture and a reduced misfit for the distributed slip model (with an elliptical shape and tapering to the fault plane edges for the seismic slip) for the 2009 L'Aquila earthquake compared to the uniform slip model on a rectangular fault plane. Thus, source parametrizations with a tapering to the edges such as elliptical source models can be a more favourable shape for a very shallow dislocation source. As we have used only synthetic data, this problem is out of the scope of this study. Nevertheless, this issue should be solved in future work. For example, we plan to use the approach of Zhu & Rivera (2002), where the sources consist of distributed point sources along the fault plane; in this case, we could force a tapering of the point source moment at the fault plane edges to get a more realistic source model for very shallow faults.
More generally, one needs to bear in mind that the modelling of regional waveforms may be difficult in complex 3-D media with strong crustal heterogeneity. Nevertheless, Frietsch et al. (2018) showed that local/regional data can be modelled adequately with the same filter range used in this study, even in a triple junction area with prevalent volcanism and possibly more complex crustal structure and inaccuracies in the Moho depth. Moreover, other earthquake source studies have also successfully used regional waveforms with similar wave periods and at similar epicentral distances (e.g. Herrmann et al. 2011;Chevrot et al. 2011;Hayes et al. 2013;Hicks & Rietbrock 2015).
Another issue is that the weighting of the different data sets in the joint inversion is a complex problem despite efforts in the earthquake source modelling community (e.g. Sudhaus & Jónsson 2009;Funning et al. 2014;Minson et al. 2013Minson et al. , 2014b. Techniques such as multi-objective optimization could help to explore the parameter space a posteriori without any a priori weighting (e.g. Kozlovskaya et al. 2007;Chiandussi et al. 2012;Schnaidt et al. 2018). Finally, another important question in multiple fault modelling is the criterion when to stop adding more faults for the given earthquake. The use of trans-dimensional inversions with a reversible jump Markov chain Monte Carlo method (e.g. Green 1995;Sambridge et al. 2006;Bodin et al. 2012) could be helpful as the number of faults would be itself an unknown determined by the data.

C O N C L U S I O N
This study presented a new multiple fault source inversion scheme using five different data types (teleseismic P waves, teleseismic S waves, teleseismic surface waves, regional waveforms and In-SAR line of sight displacements). Synthetic inversion tests highlighted that these data types are highly complementary; the seismic data have a good temporal resolution, which helps to constrain the time-shift in the source inversions, whereas the InSAR data help to get an accurate source location, given their high spatial resolution. The fault strike, dip and rake can be estimated in a more robust way in the joint multiple fault inversions compared to single data source inversions. Trade-offs between the fault dimensions, fault depth and slip can also be significantly reduced using the joint data modelling. Only the CLVD component is difficult to model with noise-perturbed data sets.
Synthetic inversion tests based on two-fault input models show that, as expected, the lower magnitude subevents can yield the largest uncertainties, notably for the event with downdip listricity. Some inversions based on one single data type (e.g. InSAR and teleseismic data) cannot distinguish the two faults, which clearly highlights the advantages of jointly inverting different data sets. On the other hand, regional seismic data are key to discriminate the two subfault solutions for the source model with downdip listricity. This suggests that downdip segmentation can be difficult to constrain in the absence of local/regional seismic data.
This study's synthetic tests also showed that comparing the iterative modelling of multiple faults to the simultaneous inversion for all source parameters, the iterative inversion scheme can lead to substantial errors in the source model with downdip listricity. Thus, this study discourages the use of such approach in studies of real earthquakes.

A C K N O W L E D G E M E N T S
We gratefully acknowledge support from the COST Action ES1401-TIDES, supported by COST (European Cooperation in Science and 972 M. Frietsch et al. Technology) and from the NERC project NE/N011791/1. We also thank support from a UCL Impact studentship and MF thanks support from a Royal Astronomical Society grant. We thank Nader Shakibay Senobari, Yehuda Ben-Zion, Phil Meredith, Zacharie Duputel and Piero Poli for fruitful discussions. We thank Martin Mai and two anonymous reviewers whose suggestions improved this manuscript. Zhu, L. & Rivera, L.A., 2002. A note on the dynamic and static displacements from a point source in multilayered media, J. geophys. Int., 148(3), 619-627.

S U P P O RT I N G I N F O R M AT I O N
Supplementary data are available at GJ I online.
Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper. Table S1 Comparison of the fault parameters for iterative and simultaneous two-fault inversions for the input two-fault model with along-strike branching faults. This is a synthetic test corresponding to the two-fault earthquake scenario shown in Fig. 2(b), with noise added for the four seismic data sets (teleseismic P waves, teleseismic S waves, teleseismic surface waves and regional waveforms) and the InSAR data sets. The fault parameters are the centroid longitude and latitude (for UTM zone 11T), the centroid depth, the fault length, width, slip, the fault strike, dip, rake, the double couple percentage, the time-shift to the GCMT time, the seismic moment M 0 (Nm), the magnitude, and the corresponding beach-ball (Stereo). The colours of the beach-balls correspond to the larger magnitude subevent (blue) and the lower magnitude subevent (red). Figure S1 Differences between the input and output fault parameters from synthetic simultaneous inversion tests whereby the time difference between the two subevents is successively varied from 0 to 100 s for the along-strike branching event (see the corresponding fault configuration in Fig. 2b). In this synthetic test, the two fault sources have an input time-shifted between 0 and 100 s. The differences are shown for fault strike (a), fault dip (b), fault rake (c), the CLVD component (d), seismic moment (e) and the interevent time observed (f). The colours correspond to the larger magnitude subevent (blue) and the lower magnitude subevent (red). Figure S2 Histograms of the misfit for InSAR data (a) and signalto-noise distribution for seismic data (b) for the N = 100 noiseperturbed data sets for the listric fault synthetic test. The L 2 -norm misfit is calculated between the noise-free data and the N = 100 noise-perturbed synthetic LOS displacements. Ascending track (a, left-hand panel) and descending track (a, right-hand panel). The histograms of the signal-to-noise ratio in the N = 100 noise-perturbed data sets are shown for teleseismic surface waves (b, top left-hand panel), teleseismic P waves (b, top right-hand panel), teleseismic S waves (b, bottom left-hand panel) and regional waveforms (b, bottom right-hand panel). Figure S3 Histograms of the misfit for InSAR data (a) and signalto-noise distribution for seismic data (b) for the N = 100 noiseperturbed data sets for the along strike branching synthetic test. The L 2 -norm misfit is calculated between the noise-free data and the N = 100 noise-perturbed synthetic LOS displacements. Ascending track (a, left-hand panel) and descending track (a, right-hand panel). The histograms of the signal-to-noise ratio in the N = 100 noiseperturbed data sets are shown for teleseismic surface waves (b, top left-hand panel), teleseismic P waves (b, top right-hand panel), teleseismic S waves (b, bottom left-hand panel) and regional waveforms (b, bottom right-hand panel). Figure S4 Regional waveforms used in the listric fault synthetic inversion tests without noise (red) and with noise added (blue). Each subplot shows the station name in the top left corner, the azimuth in the top right, the contribution to the overall misfit in the bottom left and the epicentral distance in degrees in the bottom right. A Butterworth bandpass filter in the range from 16.6 to 33.3 s is used. Figure S5 Surface waveforms used in the listric fault synthetic inversion tests without noise (red) and with noise added (blue). Each subplot shows the station name in the top left corner, the azimuth in the top right, the contribution to the overall misfit in the bottom left and the epicentral distance in degrees in the bottom right. A Butterworth bandpass filter in the range from 125 to 180 s is used. Figure S6 Teleseismic body waves used in the listric fault synthetic inversion tests without noise (red) and with noise added (blue). Top panel: P waves. Bottom panel: S waves. Each subplot shows the station name in the top left-hand corner, the azimuth in the top right, the contribution to the overall misfit in the bottom left-hand panel and the distance in degrees in the bottom right-hand panel. A Butterworth bandpass filter in the period range from 25 to 60 s is used for P waves and a period range from 25 to 100 s for S waves. Figure S7 Regional waveforms used in the along strike branching synthetic inversion tests without noise (red) and with noise added (blue). Each subplot shows the station name in the top left-hand panel corner, the azimuth in the top right-hand panel, the contribution to the overall misfit in the bottom left and the epicentral distance in degrees in the bottom right-hand panel. A Butterworth bandpass filter in the range from 16.6 to 33.3 s is used. Figure S8 Surface waveforms used in the along strike branching synthetic inversion tests without noise (red) and with noise added (blue). Each subplot shows the station name in the top left-hand corner, the azimuth in the top right-hand panel, the contribution to the overall misfit in the bottom left-hand panel and the epicentral distance in degrees in the bottom right-hand panel. A Butterworth bandpass filter in the range from 125 to 180 s is used. Figure S9 Teleseismic body waves used in the along strike branching synthetic inversion tests without noise (red) and with noise added (blue). Top panel: P waves. Bottom panel: S waves. Each subplot shows the station name in the top left-hand corner, the azimuth in the top right-hand panel, the contribution to the overall misfit in the bottom left-hand panel and the distance in degrees in the bottom right-hand panel. A Butterworth bandpass filter in the period range from 25 to 60 s is used for P waves and a period range from 25 to 100 s for S waves. Figure S10 Trade-off plot for 100 listric fault synthetic test inversions using only seismic data with different levels of noise added. The circles mark the respective best-fitting model. The yellow star marks the inversion's solution with the lower misfit. The light blue triangles are the input parameters for the subevent with the larger magnitude whereas the pink triangles are the input parameters for the subevent with the lower magnitude (see Table S1); the same is true for the dashed lines in the histograms. We plot the following parameters: the fault strike, dip, rake, the seismic moment M 0 , the time-shift and the compensated-linear-vector-dipole (CLVD) component. The beach-ball of the input model is shown for each subevent. Figure S11 Trade-off plot for 100 listric fault synthetic test inversions using only InSAR data with different levels of noise added. The circles mark the respective best-fitting model. The yellow star marks the inversion's solution with the lower misfit.The light blue triangles are the input parameters for the subevent with the larger magnitude whereas the pink triangles are the input parameters for the subevent with the lower magnitude (see Table 2); the same is true for the dashed lines in the histograms. We plot the following parameters: the fault strike, dip, rake, slip, the centroid longitude and latitude (for UTM zone 11T), the fault length, width, the centroid depth of the fault. The beach-ball of the input model is shown for each subevent. Figure S12 Trade-off plot for 100 along strike branching synthetic test inversions using only seismic data with different levels of noise added. The circles mark the respective best-fitting model. The yellow star marks the inversion's solution with the lower misfit. The light blue triangles are the input parameters for the subevent with the larger magnitude whereas the pink triangles are the input parameters for the subevent with the lower magnitude (see Table S1); the same is true for the dashed lines in the histograms. We plot the following parameters: the fault strike, dip, rake, the seismic moment M 0 , the time-shift and the compensated-linear-vector-dipole (CLVD) component. The beach-ball of the input model is shown for each subevent. Figure S13 Trade-off plot for 100 along strike branching synthetic test inversions using only InSAR data with different levels of noise added. The circles mark the respective best-fitting model. The yellow star marks the inversion's solution with the lower misfit. The light blue triangles are the input parameters for the subevent with the larger magnitude whereas the pink triangles are the input parameters for the subevent with the lower magnitude (see Table S1); the same is true for the dashed lines in the histograms. We plot the following parameters: the fault strike, dip, rake, slip, the centroid longitude and latitude (for UTM zone 11T), the fault length, width and the centroid depth of the fault. The beach-ball of the input model is shown for each subevent. Figure S14 Trade-off plot for 100 along strike branching synthetic test inversions using InSAR data and seismic data with different levels of noise added. The circles mark the respective best-fitting model. The yellow star marks the inversion's solution with the lower misfit. The light blue triangles are the input parameters for the subevent with the larger magnitude whereas the pink triangles are the input parameters for the subevent with the lower magnitude (see Table S1); the same is true for the dashed lines in the histograms. We plot the following parameters: the fault strike, dip, rake, slip, the centroid longitude and latitude (for UTM zone 11T), the fault length, width, the centroid depth of the fault, the time-shift and the CLVD component. The beach-ball of the input model is shown for each subevent.