A Stochastic source model for the 2015 Mw 7 . 9 Gorkha , Nepal , Earthquake using Multi-Dimensional Ensemble Empirical Mode Decomposition technique

The present study aims at developing a new strategy to model the spatial variability of slip on the rupture plane using multi-dimensional ensemble empirical mode decomposition (MEEMD) technique. Here, the earthquake slip distribution is split into finite number of empirical modes of oscillation called the intrinsic mode functions (IMFs). This help in identifying the fluctuation component and trend in the slip data. The trend is positive and characterizes the nonstationary mean of the slip distribution. The fluctuation component can be modelled as a stationary random field using an exponential power spectral density function. The trend can be modeled as an elliptic patch. This new technique is demonstrated for the slip distribution of the recent Nepal Earthquake, 2015. It is observed that the new model can be used to simulate the spatial complexity of slip distribution of any earthquake.


Introduction
Finite-fault inversion provides a key to understand the earthquake source complexities which is essential to predict realistic ground motion for a scenario earthquake.The derived finite-fault rupture models from inversion contain information on fluctuations of slip on the rupture plane.In addition, these earthquake slip distributions are inherently nonstationary due to the fact that slip vanishes towards the buried edges of the fault.This clearly states the level of complexity which cannot be modelled by simple deterministic functions.In such situations, stochastic approaches provide an alternate for characterizing the spatial complexity of slip on the rupture plane.Research has been carried along this direction, over the last two decade (Somerville et al., Mai and Beroza, Lavallée et al., Raghukanth and Iyengar, Raghukanth, Raghukanth and Sangeetha [1][2][3][4][5][6]).The method for investigating earthquake slip distribution in these works are based on stationary random processes models with Gaussian properties.Limitation with these studies are, the developed models are based on Fourier analysis.But, it is well known that the PSD exists if and only if the signal is a wide-sense stationary process.For nonstationary process PSD does not exist.In such situations, Fourier-based spectral analysis which uses sine and cosine functions as the basis, in analyzing nonlinear and nonstationary data can yield distorted and incomplete information on the nature of the nonstationary slip distribution.
Huang et al. [7] pointed out that the only way to represent the physics of nonlinear and nonstationary processes is through adaptive basis functions.They have demonstrated the superiority of Empirical mode decomposition technique (EMD) technique to others by analysing wind speed, tidal velocity and earthquake data.The technique has been previously applied to understand strong motion data of past earthquakes (Huang et al.,Loh et al.,Zhang et al.,).Mallikarjuna et al. [12] has used the EMD approach to forecast the air traffic by modelling the fluctuation part using artificial neural network (ANN) technique and the trend through simple regression analysis.Raghukanth [5] has proposed a new data-driven technique, namely Bi-dimensional empirical mode decomposition (BEMD) to characterize earthquake slip distribution of past earthquakes.Here the slip distribution is modelled as a sum of fluctuation and trend component.The fluctuation component is simulated as an exponential random field.The shortcoming with this work is that modelling of the trend component has not been explored in detail.The trend component is deterministic for a particular earthquake, hence it has to be modelled as 2D nonstationary function and cannot be taken as mean value of trend part as considered in Raghukanth [5].In this study, an attempt has been made to model the trend of the slip data.[13][14][15][16][17][18][19]) using seismological (teleseismic and strong motion) and geodetic (high-rate GPS, static GPS, InSAR,ALOS-2) observations of both the mainshock and its largest aftershock.Nevertheless, we restrict our study to the modelling and analysis of finite fault slip model proposed by Yagi and Okuwaki [14] (see Figure 1).The source model indicates a north dipping low-angled thrust fault in MCT with strike/dip/rake angle of 285/10/58.6.It can also be observed from the figure that the rupture extends about 160km along the strike and 88km along the dip direction with an average rupture velocity of 3km/s.A maximum slip of about 7.5m is concentrated at a radial distance of 50km to the east of epicentre (28.147 °N lat, 84.708°E long).The total seismic moment is 9.1 × 10 20 Nm corresponding to Mw 7.9.The focal depth and depth to the top of fault are 15km and 4.58km respectively.The hypocentre of the event is located at a distance of 140 and 60km along the strike and dip directions, respectively.Rise time, which is defined as the time taken by each subfault to rupture completely is 0.8s while the source duration is 50s.The entire rectangular fault in divided into 220 subfaults of size 8km x 8km.

MEEMD
The present study introduces an advanced adaptive data analysis technique, namely the multi-dimensional ensemble empirical mode decomposition (Wu et al., [20]) technique, to analysis the spatial variability of the earthquake slips.This MEEMD method has been formulated based on the EMD (Huang [7]) and its more robust version EEMD (Huang and Wu, Wu and Huang [21][22]).In the EMD approach 1-D data, X(t) is decomposed into finite (n) modes of oscillation with certain frequencies called the intrinsic mode functions, IMFs such that where Rn (t) represents the residue of data which can either be a trend (monotonic function with one maxima) or a constant value.The EMD/EEMD approach finds its application in the recently developed MEEMD method.The theory behind MEEMD is to visualize a multidimensional dataset as a combination of data in each dimension, which is followed by the application of EEMD to slices of data in each and every dimension and finally reconstructed by combining components of similar scale along each dimension.

Fig.2. Schematics of MEEMD
A schematic flowchart explaining the MEEMD method for a two-dimensional spatial data such as earthquake slip distribution f (x, y) is shown in figure 2.
To start with, we assume the 2-D slip data of Nepal earthquake, f (x, y) as a set of 1-D slices in x-direction (i.e.along the strike).Now, each 1-D slice is decomposed using EEMD resulting in a 2-D pseudo IMF like component g(x, y).This new 2-D data is further decomposed, but this time, the data is considered as a set of 1-D slices in y-direction (i.e. down dip direction) at each location along the x-direction.These 1-D slices of g(x, y) upon EEMD decomposition results in collection of 2-D components of similar scales namely h(x, y).The derived 2-D components are further combined into a reduced set of final IMF components based on the minimal scale combination principle.

Intrinsic Mode Function
The slip distribution of Nepal earthquake has been decomposed into three dynamic IMFs each evolving around a specific frequency.With the help of these extracted IMFs, one can reconstruct the original slip data as R(x, y) in equation ( 2) refers to the residue of the slip data which is similar to R(t) in equation (1). Figure 3 depicts the extracted IMFs along with the residue and slip data.A closer look at the IMFs in figure shows that the number of extremas decreases with increasing hierarchy of IMFs.The statistics of these IMFs are presented in Table 1.The 5x5 correlation matrix clearly shows that there exists a strong correlation between data and IMFs.From figure 3, it can also be observed that the residue of the data is positive and slowly varies about the long term average.One may interpret this as a nonstationary component about which the spatial variability of the slip oscillates.This residual component which represents trend of the slip data (as in equation 3) has the highest percentage variance of 74.5% and exhibits the highest correlation of 0.9.On the other hand the fluctuation component of slip which is defined as the summation of all IMFs (as in equation 4) exhibits a percent variance of 25.5%. Trend: ( , )= ( , ) T D x y R x y

Modelling the trend component
The trend represents a positive mean value of the slip field.This nonstationary characteristic of the trend can be modelled by fitting a simple deterministic function.In this study a single elliptic patch model is used to model the trend component.This method has been previously used by Vallée and Bouchon [23] to model entire slip data.In the present, we apply the same technique but, instead of modelling the entire slip data, we try to model only the trend component.

Single elliptic patch model
As a first step, we estimate the effective rupture area of the trend based on 5%-95% Husid plot which turns out to be 128 x 68 sq.km.This effective rupture patch is shown in figure 7 as a black rectangle.Now, we model this effective rupture area as a single elliptic patch with a constant rise time constant, dk and rupture velocity, Vr.
The precise parameterization of trend component of slip is detailed in Figure 8 The final slip field can be obtained as a weighted sum of the fluctuation and trend components.In this case, explained variance of each component is taken as weightage.Figure 10 shows the final simulated slip field for the 25 th April 2015 Nepal earthquake.The simulated slip field is in good agreement with original slip data shown in figure 1.The total moment of the simulated slip is comparable to total moment of the data.random field from the spectral representation method of Shinozuka and Deodatis [24] or Fourier integral method of Pardo-Igu´zquiza and Chica-Olmo [25] using the estimated correlation lengths 8. Simulate tend component as a single patch rupture model using the parameters a, b, α, Dp for a given event.9. Obtain the final random field as a weighted sum of fluctuation and trend component.10.Fix hypocentre of the event.
An ensemble of slip fields can be simulated using the proposed methodology.The simulated slip field from this approach finds its application in the synthesis of several realizations of accelerograms using any of the analytical approaches.

Discussion and conclusion
A novel approach to model spatial complexity of slip on the fault plane using MEEMD technique has been explored in this article.The decomposition of the fault slip into IMFs brings out certain interesting features of the slip distribution of Nepal earthquake.The most interesting feature that has been observed is that, MEEMD sorts the spatial frequency content of the slip distribution into a finite set of IMFs with the highest and lowest frequency corresponding to the first IMF and trend, respectively.One can try modelling the randomness associated with each IMF separately, however, this may lead to an increase in the number of model parameters whereby reducing the efficiency of the developed model.To overcome this limitation, in the present study, it is proposed that the slip field can be modelled as a combination of stationary random field (fluctuation) and a nonstationary trend.Given a particular fault and magnitude, the proposed model can handle uncertainties in the slip distribution.The developed stochastic source model together with seismological method or Spectral Finite Element Method, can be used to synthesize an ensemble of ground motion time histories for any future earthquake and for regions lacking strong motion data.Further, randomness in the simulated time histories can be attributed to the fluctuation component while peak ground displacement and peak ground residual displacement being associated with the trend term.Another advantage with the developed model is that the estimated slip field from this approach finds its real time application in generation of seismic or tsunami hazard maps.

4 )
This statistical measure indicates the relative contribution of each component (fluctuation and trend) to the data.It is clearly observed that the fluctuation component characterizes the small scale variability in the data.Contour maps of both these components for the Nepal earthquake are shown in figure 4.

Fig. 4 .
Fig.4.Fluctuation and trend of the slip field shown in figure 1.The fluctuation component is oscillatory and can be modelled as a random filed.In the second-order analysis, variation of random field at two different locations is characterized either in space by an auto-correlation function (ACF) or in the wave number domain by a power spectral density (PSD).Power spectral density of the slip field and fluctuation are shown in figure 5. Here, the wavenumber ranges between -π/dx to π/dx.The spatial frequency content of the fluctuation component is found to be 0-0.2km −1

Fig. 5 . 2 (
Fig.5.PSD of the slip field and fluctuation component The article presents a new strategy wherein, one can attempt modelling the slip distribution as a function of fluctuation and trend component instead of modelling the entire slip which is highly erratic.

Fig. 6 .
Fig. 6.PSD at the cross sections ky = 0 and kx = 0 for the fluctuation component of Nepal event shown in figure 5.

Source model of Nepal earthquake, 2015 The 25 th
April 2015 (Mw 7.9) earthquake triggered many broadband instruments across the globe.Researchers have reported this event as the largest to have occurred in the Nepal region in the last eight decade.A multitude of inversion studies have been carried out (Hayes et al., Yagi and Okuwaki, Lay et al., Tung and Masterlark, Yue et al., Zhang et al., McNamara et al., Table1.Correlation matrix of IMFs of slip data.