1 Introduction

At 08:02 on April 20, 2013 (Beijing time), a M S7.0 earthquake (the moment magnitude calculated by the USGS was M W6.6) took place in Lushan County, near the city of Yaan in Sichuan Province, China. This earthquake is a typical blind reverse-fault earthquake and the Shuangshi-Dachuan fault of the southern segment of the Longmenshan thrust belt may be its seismogenic fault (Xu et al. 2013). Many seismologists have done much works concerning its mechanisms, e.g., source rupture process inversion using teleseismic data (Hao et al. 2013; Zhang et al. 2013b, 2014; Wang et al. 2013; Liu et al. 2013), source inversion using GPS data (Jiang et al. 2014), the analysis of GPS and seismic data near the fault (Wu et al. 2013; Jiang et al. 2013; Du et al. 2013), seismic moment tensor inversion (Lin et al. 2013), the relocation of aftershocks (Han et al. 2014; Su et al. 2013; Zhao et al. 2013; Fang et al. 2013; Zhang et al. 2013c), and strong ground motion analysis associated with acceleration record (Xie et al. 2014; Wen and Ren 2014; Mooney and Wang 2014).

Usually, spatial distribution of coseismic displacements denotes the deformation of ground surface due to an earthquake. It represents a contribution of an earthquake to the change of landforms, which may be due to local non-elastic effects or elastic deformation caused by a fault. Accurate estimates of the coseismic displacement are a worthy goal. Nowadays, there are many methods being applied to such research. The use of GPS observations is the most common and direct method (Gu et al. 2009; Jiang et al. 2013; Wu et al. 2013). However, estimates made with this method are affected by the limited number of observational stations. Coseismic displacements deduced from GPS observation are scattered and sparse, thus the details of the ground motion near the fault are not well resolved. Near-fault strong ground motion records provide another method for examining this coseismic displacement (Jin and Wang 2013; Jin et al. 2014; Chen and Loh 2007; Hu et al. 2007; Rupakhety et al. 2010; Jafarzadeh et al. 2009). In addition to the disadvantage of limited observation stations, such estimates are also affected by the correction method required by strong motion records; expert judgment is needed to assess the reasonability of each estimate (Peng et al. 2011). Satellite remote sensing images can also be used for coseismic displacement estimates (Zhang et al. 2007; Xu et al. 2014; Wan et al. 2008) though the resolution is limited. The reasonable simulation of the displacement field near the fault can provide a detailed distribution of coseismic displacement. Not only it provides overall distribution information, but also its resolution is high as well. Therefore, this study prefers to use the simulation of the wave field for estimating the distribution of coseismic displacement.

A quasi-static elastic wave equation, which omits the acceleration, is often applied for simulations of coseismic displacement. In fact, coseismic displacement is not a kind of independent displacement field, rather it is the final result of the development of the ground motion displacement. Therefore, I suggest that a simulation of dynamic displacement can provide more accurate information on coseismic displacement. This paper uses the elastic wave equation directly and applies the spectral element numerical method to simulate the displacement field of the Lushan earthquake due to source rupture process. The coseismic displacement is estimated by the permanent value of the displacements. Prior to the displacement simulation, the simulated ground velocities were compared with the velocities from strong ground motion records. Their agreements confirm the reasonability of the source and media model. Furthermore, the comparisons between the simulated results and coseismic displacements deduced from GPS validate our estimates. In Sect. 5, the distributions of coseismic displacement are analyzed and discussed, highlighting some interesting and significant phenomena concerning the ground motion observed in this region.

2 Simulation method

The spectral element method (Komatitsch and Vilotte 1998; Komititsch and Tromp 2002; Komititsch et al. 2004) was applied in the simulation. The underground propagation in a given medium is described by the wave equation with boundary and initial conditions that can be written as

$$\begin{aligned} \rho {\user2{\ddot{u}}} = \nabla \cdot {\sigma} + {\user2{f}} \hfill \\ {\sigma} = {C}:{\varepsilon} \hfill \\ {\varepsilon} = \frac{1}{2}\left[ {\nabla {\user2{u}} + (\nabla {\user2{u}})^{\text{T}} } \right], \hfill \\ \end{aligned}$$
(1)

where \({\user2{u}}\) is the displacement vector, \({\sigma}\) is the stress tensor, \({C}\) is the stiffness tensor, \(\rho\) is the density, \({\user2{f}}\) is the external force, and \({\text{T}}\) denotes transposition. The colon symbol in the equation denotes multiplication by a tensor and the dots above the \({\user2{u}}\) denote a time derivative.

For the above equation, an arbitrary test function, \({\user2{w}}\), is multiplied on both sides. Using integration by parts, Eq. (2) is obtained as follows:

$$\int_{\varOmega } {\rho {w} \cdot {\user2{\ddot{u}}}{\text{d}}\varOmega + \int_{\varOmega } {\nabla {w}:{C}:\nabla {\user2{u}}{\text{d}}\varOmega = \int_{\varOmega } {{w} \cdot {\user2{f}}{\text{d}}\varOmega } } } + \int_{\varOmega } {{w} \cdot {\user2{t}}{\text{d}}\varGamma } ,$$
(2)

where \(\varOmega\) is the computation zone, \(\varGamma\) is the boundary that includes the free surface and the artificial boundary, and \({\user2{t}}\) is the traction on \(\varGamma\).

For the spectral element method (SEM), the area to be calculated is first divided into non-overlapping discrete small elements with six surfaces. The Gauss-Lobatto-Legendre (GLL) points are computed for each element. Finally, denoting the global vector of unknown displacement by \({\user2{U}}\), the matrix equation of Eq. (2) was obtained as follows:

$${\user2{M{\ddot{U}}}} + {\user2{KU = F}},$$
(3)

where \({\user2{M}}\) is the mass matrix, \({\user2{K}}\) is the stiffness matrix, and \({\user2{F}}\) is the source term. In this procedure, the application of GLL points in conjunction with the GLL integration rule renders the mass matrix \({\user2{M}}\) exactly diagonal, thereby resulting in drastic reductions in computational cost of the algorithm. This is a major advantage of the SEM. In the following simulation, the time interval is set as 0.01 s and the small elements are set as 2 km. In each direction of small elements, 7 GLL points are applied.

3 Source and media models

Figure 1 shows the region of latitude 29.5°N–31.5°N and longitude 102°E–104°E. In the figure, the study area is indicated by a yellow square having sides of 110 km. The media parameters applied in our simulation were constructed by an interpolation of the crust1.0 model with the topography. Figure 2 shows the varying Vp along the profile AB shown in Fig. 1. Obviously, the model used in computation is 3D irregular layer media. The rupture process of the Lushan earthquake has been investigated by many inversion works. On the whole, the results of these inversions are similar; however, their differences are not trivial (Hao et al. 2013; Zhang et al. 2013a, 2014; Wang et al. 2013; Liu et al. 2013; Jiang et al. 2014). Zhang et al. (2013b, 2014) provided not only the slip distribution but also the slip-rate time function on the faults. I simulated ground motion using the sources inverted by Zhang et al. (2013b) and Zhang et al. (2014) and compared them with observed records respectively and found those of Zhang et al. (2013b) providing a closer agreement for the media model used in this study. So their result (Zhang et al. 2013b) was applied as the source function in this study. The details of the source can be found in Zhang et al. (2013b). The left panel of Fig. 3 shows the slip distribution on the fault and the right panel displays the slip-rate time function along the strike across epicenter.

Fig. 1
figure 1

Simulation region and observed stations

Full size image
Fig. 2
figure 2

Velocity structure of P wave under the profile AB

Full size image
Fig. 3
figure 3

Left panel the earthquake moment distribution on the fault surface. Right panel the slip rate along strike across the epicenter

Full size image

4 Comparisons of simulated results with strong ground motion records

In Fig. 1, near the fault there are three observed stations of strong ground motion denoted by red open circles and six GPS stations denoted by red triangles. The epicenter is denoted by the red star and the red box is the projection of the fault plane on the ground surface. The left column in Fig. 4 shows the acceleration records in three directions at stations 51BXY, 51BXM, and 51LSF. Obviously, they contain abundant high-frequency information. Due to the low frequency of the source (less than 0.1 Hz) and the limitation of computation technique and computation equipment, only the acceleration in low-frequency band can be simulated. Thus, the low-frequency simulation results and acceleration records need to be filtered by the same low-pass filter for comparison. The middle column in Fig. 4 shows the comparison between the simulation results and the acceleration records filtered by the same low-pass filter. The right column shows the comparison between the simulation velocities and the recorded velocities which were also filtered by the same low-pass filter. Comparisons of the 18 records at three stations show that both velocity and acceleration simulation results have waveforms and values similar to those of the observed records. For instance, both have a large trough and peak with similar periods. These agreements illustrate that the source used can represent the actual main rupture of the M S7.0 earthquake and that the wave propagation in the model can reflect the main wave information at low frequencies. Thus, the estimation of coseismic displacement with lower frequencies based on the same source and media model can be considered trustworthy.

Fig. 4
figure 4

Acceleration record in three directions at 51BXY, 51BXM, and 51LSF ground motion stations, and the comparison of velocity records and acceleration records with that of simulated results

Full size image

Usually, the displacement has a lower frequency than the velocity. Figure 5 shows the simulated displacements at stations 51BXY, 51BXM, and 51LSF. The unfiltered result is denoted by black lines and the filtered one is represented by a red line. Obviously, the displacements in all three directions have static displacements that can be considered as the actual permanent movement at the stations. To remove probable high frequencies in the static displacement, we used the filtered value to estimate the coseismic displacement field of the Lushan earthquake.

Fig. 5
figure 5

Comparisons of the filtered and unfiltered simulated displacements at stations 51BXY, 51BXM, and 51LSF

Full size image

Before analyzing the coseismic displacement field, a further evaluation of the simulated coseismic displacement needs to be undertaken by comparing GPS stations. We gathered GPS data from six stations within the computation domain. Their horizontal and vertical displacements, as compared with the simulated displacement at the corresponding locations, are shown in Table 1. The absolute difference (about 2 cm) of point LS05 is the largest. However, the difference is not large relative to GPS measurements of the point (8.75 cm). The simulated displacements at LS07 are mostly consistent with displacement records; their differences are less than the observed GPS errors. The difference in other measurements at points LS01, SCTQ, and BXB is less than 1 cm. The relative difference at station QLAI is the greatest, possibly because of its special position. The following Sect. 5 will show that this station is located very close to the fault line AB, where very strong changes in ground motion occurred. Despite this, the relative difference at station QLAI is nearly equal to its relative measurement error of GPS. Therefore, the differences between the simulation results and GPS data are considered reasonable, and the simulation results can represent the actual coseismic displacement. Figure 6 shows the comparisons of horizontal displacement vectors at stations LS01, LS07, and BXB. Not only are their values similar, but also their directions of motion agree to each other. The vector comparisons on the other three stations are not as good as those on the three stations LS01, LS07, and BXB, so their comparisons are omitted. In Figs. 7, 8, 9, and 10, there is a starting (*) point at 30.20°N, 102.94°E, which has an upward moving displacement of 19.843 cm from level measurement (Hao et al. 2014), whereas the corresponding simulation result is about 20 cm. Their difference is, thus, less than 2 mm. The comparisons above prove that our estimates of coseismic displacement indeed represent actual ground motions.

Table 1 Comparisons between the simulated results and coseismic displacement determined by GPS
Full size table
Fig. 6
figure 6

Comparison of horizontal vector motions at the three GPS stations and the simulated results

Full size image
Fig. 7
figure 7

Distribution of vertical coseismic displacements. Asterisk indicates the level measurement point. The blue box denotes the projected zone of the fault surface on the ground, the triangles denote station

Full size image
Fig. 8
figure 8

Distribution of north-south-oriented coseismic displacement. The blue box denotes the projected zone of the fault surface on the ground, the triangles denote station

Full size image
Fig. 9
figure 9

Distribution of east-west-oriented coseismic displacement. The blue box denotes the projected zone of the fault surface on the ground, the triangles denote station

Full size image
Fig. 10
figure 10

Distribution of horizontal vectors of coseismic displacement. The blue box denotes the projected zone of the fault surface on the ground

Full size image

5 Coseismic displacement distribution analysis

Figures 7, 8, 9, and 10 show the distribution of the coseismic displacements obtained from the static displacement of simulated wave propagation. To conveniently describe displacements, we introduce two concepts: one is the projection region that corresponds to the projected zone of the fault surface on the ground as denoted by a blue box in Figs. 7, 8, 9, and 10; the other is the fault line AB (Figs. 7, 8, 9, 10) corresponding to the fault trace on the surface. The star symbol indicates the epicenter.

Figure 7 illustrates the distribution of vertical displacement. An uplift is obvious in the vicinity of the epicenter, where upward motion is more than 25 cm. The uplift spreads in a circular pattern as indicated by 5-cm interval contour lines with a radius <25 km. The appearance of the uplift agrees to the assumption of Chen et al. (2014) who showed that there should be an uplift zone on the hanging wall. In the northwest, away from the projection region, there is a zone of depression with a decrease in elevation of <1 cm. This region is produced by the uplift motion. Along the fault line AB, there are small areas distributed like beads with large upward motion that exceeds 20 cm, and near the south end of AB, there is also a weak depressed zone with upward motion around it. Repeatedly, near the epicenter in this figure there is a star point (30.2°N, 102.94°E) with an upward motion of 19.843 cm (Hao et al. 2014), matching our simulated value of 20 cm. This similarity shows the reliability of our result. In comparison with the hanging wall, the motion on the footwall is simple with a rising displacement less than 1 cm. On the whole, the uplift near the epicenter is dominant that controls the vertical motion near the fault.

The distribution of motions in the north-south direction is shown in Fig. 8. It is clear that the movement on the hanging wall can be divided into three parts: outside the projection area, southward motion dominates with displacement less than 3 cm; within the projection area, motion in the southern part is toward the south with displacement more than 7 cm, while in the northern part it is toward the north and less than 2 cm. These opposite motions result from uplift near the epicenter. The rising ground pushes the material adjacent to the epicenter to the north and south. Note that along the fault line AB motions are very complex. The northward and southward motions appear alternatively. This alternating motion is apparently determined by the alternative rise and fall along the fault line AB. In the northern section of the line, the maximum motion to the north exceeds 8 cm. The northward motion near the left boundary of the simulation zone may be explained by the downward-oriented motion near this area. Of course, the amplitude of this motion is small and is only a subordinate motion. On the footwall, the northward-oriented motion only occurs in a half circle near the fault with a value less than 2 cm.

Figure 9 shows the contours of east-west-oriented motion. On the whole, motion on the hanging wall is >0 cm, indicating move to the east, whereas on the footwall the value is <0 cm, denoting move to the west. The distribution also shows that in the fault projection area the movement is complex, as demonstrated by a strong variation of the contours. In the northeastern part of the projection, there are two strong movement zones: one lies just to the east of the epicenter with the largest motion exceeding 14 cm, and the other appears near the northern section of fault line AB and exceeds 16 cm. Near the center of the projection, there is a small zone where the motion in the west-east direction is close to zero; this region lies near the epicenter and has an area of about 3.0 km × 3.0 km. In contrast to the hanging wall, the footwall shows not only low level of motion (with magnitudes <3 cm), but also a weak change of motion.

To illustrate the horizontal motion more clearly, vector horizontal displacements are shown in Fig. 10. This illustrates very interesting movements. On the hanging wall, the epicenter and fault line AB control the distribution of ground motion. Outside the projection area, the ground moves toward the projection. In the center of the projection, the epicenter controls the motion. Although the epicenter itself has little horizontal displacement, it forces the surrounding zone to move toward the northeast, east, and southeast as the epicentral distances increase. When these motions intersect the fault line, they stop or change direction. Therefore, in the northern section of fault line AB, the motion becomes strike slip, whereas in the middle section it is perpendicular to AB and in the southern section it is oblique to AB. Differently on the footwall, within 40 km from the fault line, the slight movements are toward the northwest. Furthermore, in the northeast and southwest parts of the fault projection, there appear to be some rotational motions, such as eastward motions turning toward southeastward and then southward. On the whole, the southeastward motion is the main horizontal movement observed. From the above analysis, we can see that southeastward and upward motion dominates in the area, whereas the northward movement is driven by upward motions.

6 Conclusions

The paper simulated the wave field near the fault that hosted the M S7.0 Lushan earthquake using the spectral element method. The inverted rupture process has been applied as the source. The media model was constructed based on the parameters of the crust1.0 model. First, the simulated ground velocity was compared with the ground motion records that were filtered to a low-frequency band. The agreement in this comparison proves that the media and source applied in the computation can provide a reasonable simulation for the low-frequency band. Therefore, the static displacement in the simulated wave field can be used to estimate coseismic displacement. The estimated coseismic displacements were then compared with the corresponding GPS data. These comparisons indicate that the displacements in horizontal or vertical direction are similar. In Sect. 5, we carefully analyzed the distribution of simulated coseismic displacement and found that the distribution suggests the following characteristics.

  1. (1)

    The motion on the hanging wall is greater than that on the footwall.

  2. (2)

    For vertical motion, the greatest displacement (>20 cm) occurs close to the epicenter on the hanging wall, manifesting as an uplifted zone. Twenty-five kilometers away from this location, the vertical motion decreases rapidly to less than 5 cm.

  3. (3)

    On the hanging wall, east-west-oriented motion is all toward the east. The largest motion (about 16 cm) occurs in the northern section of the fault.

  4. (4)

    For south-north-oriented motion, there are two oppositely moving regions: one toward the north and the other toward the south. They are controlled by upward motion proximal to the epicenter.

  5. (5)

    The vector distribution of horizontal motion shows that the horizontal motion is complex. Rotational and along-strike movements all appear at different epicentral distances. In contrast to the hanging wall, motion on the footwall is simple with displacement <1 cm.

The above analysis shows that an estimation of coseismic displacement from a simulated wave field provides more information on ground motion for the Lushan earthquake. This is meaningful to help us understand the topography change resulting from a strong earthquake. However, the above descriptions only prove that the estimate is effective for the Lushan earthquake. Further study is needed to clarify whether this kind of estimation can be applied to other earthquakes or not. Besides, in Lushan earthquake there were lots of coseismic landslides; do the coseismic landslides affect the coseismic surface displacement? What relations are there between the coseismic landslides and the coseismic displacement? These scientific questions need to be answered in the future.