1 Introduction

Complex sedimentary basin sites have a significant amplification effect on ground motion. Research has shown that the basin effect is an important factor that causes earthquake disasters and localized anomalies of ground motion. As a typical complex site condition, sedimentary basins may aggravate earthquake disasters in many cases, such as the Mexican earthquake (Eduardo 1999), the Northridge earthquake in the USA (Graves et al. 1998), and the Kobe earthquake in Japan (Kawase 1996).

Numerical simulation, a common method to study the basin effect, can simulate the process of seismic wave propagation in a medium with any complex underground structure and calculate the ground motion. Kawase and AKI (1989) obtained similar results with the 1985 Mexico earthquake records by simulating a two-dimensional basin model. A series of three-dimensional numerical simulation works were carried out in the Los Angeles Basin (Graves (2000) and Olsen (2000)) and Santa Clara Valley, CA (Harmsen et al. (2008)), in which different source parameters and models were compared and analyzed. In addition, Wald and Graves (1998) in Kobe; Satoh et al. (2001) in Sendai; Kagawa et al. (2004) in Osaka; Shiann Jong Lee et al. (2008) in Taipei; and Frankel et al. (2009) in Seattle and Washington performed similar three-dimensional numerical simulation work.

A mature numerical simulation method is to build an integral model including the source, which can comprehensively study the influence of the source parameters, propagation path, and site structure of large earthquakes on strong ground motion. Galvez et al. (2014) simulated the source rupture process of the 2011 Japan m 9.0 Earthquake and reproduced two important observation features; Savran and Olsen (2019) simulated the Chino Hills earthquake in California in 2008 and studied the influence of the uncertainty in velocity and density in the structure model on the simulation results. Graves et al. (2011) of the Southern California Earthquake Center (SCEC) built a cybershake computing platform and promoted the application of numerical simulation methods in seismic risk analysis.

Due to the large scale of simulation models, versions including the source have low accuracy in describing the precise structure and terrain of localized sites when limited by the calculation capacity and have limited ability to simulate high-frequency ground motion. Therefore, many scholars have carried out large-scale numerical simulations of ground motion on supercomputing platforms to obtain a higher simulation accuracy and cutoff frequency. Roten et al. (2016) simulated a 7.7-magnitude earthquake south of the San Andreas fault on NCSA Blue Waters and OLCF Titan with a cutoff frequency of 4 Hz and a spatial resolution of 25 m; Fu et al. (2017) conducted a large-scale nonlinear earthquake simulation of the Tangshan earthquake on Sunway TaihuLight with a cutoff frequency of 18 Hz and a spatial resolution of 8 m.

Large-scale earthquake numerical simulation requires plentiful computing resources. To consider the influence of multiple potential sources and the uncertainty of site structure, it is necessary to carry out multiple calculations in the process of numerical simulation research or engineering seismic analysis for a specific site. To accomplish multiple simulations with a high cutoff frequency and spatial resolution with limited computing resources, a feasible method is to build a localized site model without a source and solve the scattering problem of the infinite domain. In this paper, the artificial boundary and the explicit finite element method are combined to develop a multitransmitting formula (MTF) in automatic dynamic incremental nonlinear analysis (ADINA) software; the input and scattering of seismic waves on the boundary are achieved, and the accuracy of the calculation results is verified. This method can obtain a “pure” site effect independent of the source and propagation path and meet the needs of basin effect research and major engineering input ground motion; the model is applied to the numerical simulation of the basin effect in the Yuxi Basin.

2 Method

The explicit finite element numerical simulation method has great advantages in the study of complex systems. Although it is conditionally stable, it has the advantages of low computational complexity and high efficiency in general. Xiaojun et al. (1992) proposed the method of combining the central difference with the Newmark constant average acceleration with the second-order accuracy, which improved the accuracy of calculation. Xueliang and Yifan (2002), Xueliang et al. (2006) combined Taylor expansion with the generalized Newmark-β method, gave an explicit step-by-step integration scheme with second-order accuracy, and realized the direct input of acceleration in the explicit algorithm.

The MTF is a kind of artificial boundary put forward by Liao Zhenpeng in 1984. The basic idea of the MTF is to assume an artificial wave velocity near the boundary with a direction along the normal of the boundary and then simulate the physical process of the wave transmitting outside the computational region while simulating multitime systems with the error wave method to improve accuracy. The MTF is universal because it is independent of the control differential equation of the infinite media and the physical character of the boundary; moreover, the MTF easily matches the step-by-step calculation of the finite element and finite difference methods, and the precision of the MTF is controllable (Liao 1996). The MTF is as follows:

$$ {u}_0^{p+1}=\sum \limits_{j=1}^N{\left(-1\right)}^{j+1}{C}_j^N{u}_j^{p+1-j} $$
(1)

in which:

$$ {C}_j^m=\frac{m!}{\left(m-j\right)!j!} $$
(2)

By combining the MTF with the lumped mass finite element method, a complete decoupling explicit finite element method can be developed. This method can solve precise site models with complex structures, and it has great advantages in calculation capacity and efficiency.

The numerical simulation work in this article was based on ADINA software. ADINA is a dynamic nonlinear finite element software developed by ADINA R&D Company under the leadership of Dr. K. J. Bathe; the software has advantages in pretreatment, fluid-solid coupling and nonlinear calculation areas. However, ADINA does not provide the artificial boundary that is needed in the calculation of wave propagation in infinite media, so we first incorporated the MTF into ADINA by secondary development and verified the validity of the boundary.

3 Development process

Matching with the explicit spatiotemporal decoupling MTF formula, ADINA software provides a spatiotemporal decoupling finite element algorithm utilizing the combination of the lumped mass finite element method and the central difference method. The basic motion equation is:

$$ \left(\frac{1}{\Delta {t}^2}\left[M\right]+\frac{1}{2\Delta t}\left[C\right]\right){u}^{p+1}={P}^p-\left(\left[K\right]-\frac{2}{\Delta {t}^2}\left[M\right]\right){u}^p-\left(\frac{1}{\Delta {t}^2}\left[M\right]-\frac{1}{2\Delta t}\left[C\right]\right){u}^{p-1} $$
(3)

By combining the equations of motion (1) of the boundary nodes in the MTF with the equations of motion (3) of the internal nodes in the finite element method, a decoupled numerical simulation method for the near-field ground motions can be established. The basic idea is to use the interpolation method to rewrite the MTF as a boundary node motion equation expressed by the node displacement of the finite element spatiotemporal discrete grid and calculate the response of the boundary node step by step. Therefore, the formula of the MTF formula expressed by node displacement is:

$$ {u}_0^{p+1}=\sum \limits_{j=1}^N{\left(-1\right)}^{j+1}{C}_j^N{T}_j{u}_{d,j}^{p+1-j} $$
(4)

To simulate the ground motion of the model without a source, it is necessary to input the seismic wave on the boundary. The input seismic wave generates a free field in the infinite half space of the model. Since the MTF can simulate the response of only the outward wave at the boundary point, when elucidating the motion equation of the boundary point, the free field ur must be removed from the whole wave field u to separate the outward wave us. The equation is:

$$ {u}_s=u-{u}_r $$
(5)

Let the displacement u = us in the MTF (4). The equation is:

$$ {u}_0^{p+1}=\sum \limits_{j=1}^N{\left(-1\right)}^{j+1}{C}_j^N{u}_{s,j}^{p+1-j}+{F}_0^{p+1} $$
(6)

\( {F}_0^{p+1} \) in formula (5) can be calculated from the input seismic wave. Therefore, this method can simulate the ground motion considering the input of seismic waves and the simulation of the MTF.

In this paper, we first develop a one-dimensional MTF and then extend it to two-dimensional and three-dimensional cases. In this process, the results of the one-dimensional MTF model can be used as the input of the seismic wave free field when two-dimensional and three-dimensional models are calculated. Furthermore, for the two-dimensional and three-dimensional MTF, it is necessary to establish a set of boundary node numbering methods to determine the boundary direction and relative position of all nodes. Figure 1 is the flow chart of the secondary development program of the MTF in ADINA.

Fig. 1
figure 1

Flow chart of the secondary development program of the MTF in ADINA

Full size image

4 Verification of 3D model simulation results

In this paper, the accuracy of the development results is verified by comparison with the analytical solution. The error of the calculation results of the one-dimensional and two-dimensional models is small, and the verification process is not covered here.

To verify the three-dimensional MTF, a three-dimensional cube model with a single-layer medium is built. The top of the model is a free surface, while the other five surfaces satisfy the MTF. Three components of displacement seismic waves are input at the bottom, with the direction perpendicular to the bottom of the model. Considering the size of the target area and the physical properties of the underground strata, the model size of the cube is set as 3000 m, with a density of 2000 kg/m3, a shear wave velocity of 3000 m/s, a Poisson ratio of 0.33, and a compression wave velocity of 6000 m/s. The mesh size of the model is 100 m, with 27,000 cells and 29,791 nodes. The time step is taken as 0.005 s, which is close to the time step in the actual calculation. The model is shown in Fig. 2. The model is divided into six parts: the middle area is the site structure to be calculated, and the box around and below is the MTF area.

Fig. 2
figure 2

Waveform and amplitude spectrum of the input pulse wave

Full size image

The input seismic waves of the three components are all delta pulse shear waves with widths of 1 s, and the formula is as follows:

$$ T(t)=\left\{\begin{array}{c}16{t}^3\kern9.5em 0\le \mathrm{t}<0.25\\ {}1-48\mathrm{t}{\left(t-0.5\right)}^2\kern4.25em 0.25\le \mathrm{t}<0.5\\ {}1+48\mathrm{t}{\left(t-0.5\right)}^2\kern3.75em 0.5\le \mathrm{t}<0.75\\ {}-16{\left(t-1\right)}^3\kern7.5em 0.75\le \mathrm{t}<1\\ {}0\kern11em \mathrm{t}<0\ \mathrm{or}\ \mathrm{t}\ge 1\end{array}\right. $$
(7)

For node n, the formula of its analytical solution is:

$$ {T}_N(t)=T\left({t}_1\right)+T\left({t}_2\right) $$
(8)

in which:

$$ {t}_1=t-\frac{\left(N-1\right)\ast 5}{V_s},{t}_2=t-\frac{\left(101-N\right)\ast 5+500}{V_s} $$
(9)

The waveform and amplitude spectrum of the input seismic wave are shown in Fig. 3.

Fig. 3
figure 3

Model diagram for verifying the simulation results of the 3D MTF

Full size image

The simulation results are as follows. Figure 4 shows the comparison of the displacement response between the simulation results and the analytical solution at the bottom node, and Fig. 5 shows the comparison at the free surface node. For the boundary node, the numerical results in the X and Y directions are slightly different from the analytical solutions. The main reason is the concentration of the seismic wave reflected from the free surface back to the boundary. The maximum error value is approximately 0.004 m, which is equivalent to 0.4% of the peak value of the pulse wave. The difference between the numerical results and the analytical solutions in the Z direction is small, and the maximum error is approximately 0.002 m, which is equivalent to 0.2% of the peak value of the pulse wave. A possible reason is the high P wave velocity of the model. For the free surface nodes, the numerical results of the three components are slightly different from the analytical solutions, and the maximum error value is approximately 0.004 m, which is equivalent to 0.4% of the peak value.

Fig. 4
figure 4

Comparison of numerical and analytical solutions of boundary nodes

Full size image
Fig. 5
figure 5

Comparison of numerical and analytical solutions for free surface nodes

Full size image

In general, the second development of the MTF in ADINA successfully simulates the process of seismic wave propagation from the calculation area to the infinite area. In the numerical test results, the incident wave is reflected by the free surface and then propagates out of the calculation area through the MTF. Some seismic waves are still reflected back and forth in the calculation area, and their displacement peak value is less than 0.0005 m over the time period of interest, which is less than 0.05% of the input peak value and can normally be ignored.

5 Example of Yuxi Basin and analysis

The Yuxi Basin, located in Yunnan Province, is one of the larger sedimentary basins located in the southern part of the Puduhe fault zone. It is about 10 km from east to west and 20 km from south to north. In this area, the Pliocene is composed of lacustrine, swamp, and alluvial deposits; the lower Pleistocene is alluvial lacustrine, and the middle and upper Pleistocene is dominated by alluvial proluvial and eluvial soil; the Holocene is composed of river terraces and fluvial deposits; the Quaternary is mainly alluvial proluvial gravel and clay. The east side of Pudu River fault forms two sedimentary centers, one of which is located in the north of the basin, with the maximum sedimentary thickness of about 800 m; the other is located in the middle of the basin, with the maximum sedimentary thickness of more than 800 m, which is the lowest part of the whole basin. The location and sedimentary thickness of Yuxi Basin are shown in Figs. 6 and 7.

Fig. 6
figure 6

The geographical location map of Yuxi Basin

Full size image
Fig. 7
figure 7

Distribution of sedimentary thickness of Yuxi

Full size image

Based on the buried depth model of the Yuxi Basin basement and relevant data given by He et al. (2013), data from 17 seismic exploration sections, 17 electrical exploration sections, 60 surface wave explorations, and 67 boreholes are synthesized, and a velocity structure model of the Yuxi Basin is built (Zhen et al. 2017). The depth of the model is about 800 m, which is divided into six layers, reflecting the structure of sedimentary soil layer in the basin area. Compared with the previous basin numerical simulation studies, the model structure is more refined.

The Qujiang fault near the Yuxi Basin is selected as the seismogenic fault (the Tonghai earthquake in 1970). By using the finite difference method of recurrence of velocity and stress (Guo 2016), the propagation of strong ground motion and the simulation of the fault source process are realized. The seismic wave field at the bottom of the bedrock in the Yuxi Basin is calculated as the input. The time history of input ground motion is shown in Fig. 8.

Fig. 8
figure 8

Seismic waves of bedrock at the bottom of Yuxi Basin calculated by the finite difference source model

Full size image

The peak acceleration in the east-west (E-W) direction is 1.32564 m/s2, and the maximum displacement is 0.02857 m; the peak acceleration in the north-south direction is 0.73217 m/s2, and the maximum displacement is 0.03803 m; the peak acceleration in the vertical direction is 0.55017 m/s2, and the maximum displacement is 0.03135 m, corresponding to a plane wave with the propagation direction perpendicular to the bottom surface of the model. The numerical simulation adopts displacement input. The distribution of the magnification of the displacement in the E-W, north-south, and vertical directions is calculated.

Figure 9 shows the distribution of the displacement magnification of the E-W ground motion. As shown in the figure, the amplification of E-W ground motion in the Yuxi Basin shows a certain basin focusing effect; that is, the amplification of ground motion in a specific area in the middle of the basin is usually large and sometimes small. This is due to the influence of the undulating basement of the basin and the superposition of seismic waves of different phases, which results in an amplification or reduction in regional ground motion. Among these responses, the distribution of the focusing effect that causes the increase in magnification has a certain relationship with the depth of the sedimentary layer of the basin. However, the areas with the largest magnification are not in the deepest position of the sedimentary layer. On the one hand, they are affected by the structure of the basin; on the other hand, they are due to the wave velocity and Young’s modulus of the deepest layers in the sedimentary layer of the basin, which has minor differences with bedrock and leads to an inconspicuous amplification effect. At the same time, the amplification distribution of the ground motion in the E-W direction of the Yuxi Basin also shows a certain edge effect; that is, the area located at the edge of the basin with a shallow sedimentary layer has a larger amplification than other areas. This effect is obvious in the northern edge of the basin because the depth of the sedimentary layer at the edge of the basin changes rapidly and the bedrock surface is steep; furthermore, because the sedimentary layer in this area is a shallow and low-speed layer, its wave velocity and Young’s modulus are quite different from the bedrock in the lower layer.

Fig. 9
figure 9

Distribution diagram of the amplification ratio of E-W seismic displacement

Full size image

Figure 10 shows the distribution of the displacement magnification of the north-south earthquake. As shown in the figure, the seismic amplification in the north-south direction of the Yuxi Basin shows a more obvious basin focusing effect and edge effect than in other directions. The area with an obvious focusing effect is located at the place with the deepest sedimentary layer in the middle of the basin, while the area with an obvious edge effect is located at the edge of the east and west sides of the basin, especially at the location with the steepest bedrock surface in the west of the basin, with the most obvious edge effect.

Fig. 10
figure 10

Distribution diagram of the amplification ratio of north-south seismic displacement

Full size image

In the simulation results, the displacement magnification of the velocity structure model of Yuxi Basin for the E-W earthquake is approximately 2 times, the corresponding peak displacement is approximately 0.12 m, the peak acceleration is approximately 3.8 m/s2, and the displacement magnification of most areas is between 1.2 and 1.6 times; the maximum displacement magnification for the north-south earthquake is approximately 1.9, and the corresponding peak displacement is approximately 0.11 m, the peak acceleration is approximately 2.5 m/s2, and the displacement magnification in most areas is also between 1.2 and 1.6 times. In general, the E-W amplification of the input ground motion in the Yuxi Basin is slightly larger than the north-south amplification. However, the edge effect of the north-south earthquake motion is more obvious.

Figure 11 shows the distribution of displacement magnification of vertical ground motion. As shown in the figure, the amplification effect of the basin sedimentary layer on the vertical ground motion is minimal. The displacement magnification of most areas is between 1.0 and 1.2 times, and the maximum displacement magnification is approximately 1.3 times, which is located in the deepest part of the basin sedimentary layer.

Fig. 11
figure 11

Distribution diagram of the amplification ratio of vertical seismic displacement

Full size image

Fu and Gao (2016) conducted finite difference ground motion simulation in Yuxi basin, studied the amplification effect of basin on ground motion within 2–10 s, and gave the relationship between basin amplification coefficient and equivalent sedimentary thickness. In this paper, a local site model without seismic source and a MTF boundary are used, so a more precise model is calculated by using a personal computer. The low-frequency simulation results in this paper are similar to those of Fu and Gao (2016), both of which reflect the focusing effect of basin and the influence of sedimentary thickness on seismic amplification, and the medium-high-frequency simulation results reflect the edge effect of the basin better. Figure 12 shows the distribution of basin magnification under different frequencies. As shown in the figure, the high-frequency ground motion is more affected by the local structure of the basin, while the low-frequency more affected by the overall structure.

Fig. 12
figure 12

Distribution of the amplification ratio of E-W seismic acceleration the corresponding frequency of acceleration response spectrum are 1 Hz, 0.66 Hz, 0.5 Hz, and 0.25 Hz

Full size image

6 Conclusion and discussion

In this paper, the development of the MTF in ADINA software is completed, and the calculation results are verified. The error is less than 0.5%, so this method can be used in the numerical simulation of the basin effect of a fine basin structure. The calculation example in this paper can be completed in 2 h by using a single computational node or a personal computer, which provides the basis for multiple simulation studies of different parameters of different potential sources in the same site. The cutoff frequency of this simulation is relatively low. When clusters are used and the calculation time is increased, it is expected that the calculation cutoff frequency of the small-scale model can be increased to more than 15 Hz.

The ground motion simulation of the Yuxi Basin shows a strong amplification effect on ground motion. Generally, the amplification ratio of the displacement peak in the basin is between 1.2 and 1.6 times, and the distribution is uneven, showing multiple extremum areas. The amplification ratio of the extremum area caused by the focusing effect of the basin can reach approximately 2 times.

The focusing effect of the Yuxi Basin is manifested mainly in the increase in magnification caused by the superposition of direct waves, multiple reflection waves, and surface waves and the decrease in magnification caused by the mutual cancelation of seismic waves with opposite phases in some areas. There is a certain relationship between the focus effect distribution of the E-W seismic motion and the depth of the second and third layers in the basin model, while the focus effect of the north-south seismic motion is concentrated mainly in the deepest area of the basin. This may be related to the structure of the narrow E-W and wide north-south features in the Yuxi Basin.

The simulation and recognition of the basin effect in this paper are consistent with previous research, and its distribution characteristics can provide a reference for the seismic fortification planning of related areas. We plan to use multiple potential sources and more refined basin structural models to further study the basin effect in this area.