Simulating Shallow Soil Response Using Wave Propagation Numerical Modelling in the Western Plain of Taiwan

This study used the results from 45 microtremor array measurements to construct a shallow shear wave velocity structure in the western plain of Taiwan. We constructed a complete 3D velocity model based on shallow and tomography models for our numerical simulation. There are three major subsurfaces, engineering bedrock (VS = 600 m s-1), Pliocene formation and Miocene formation, constituted in the shallow model. The constant velocity is given in each subsurface. We employed a 3D-FD (finite-differences) method to simulate seismic wave propagation in the western plain. The aim of this study was to perform a quantitative comparison of site amplifications and durations obtained from empirical data and numerical modelling in order to obtain the shallow substructure soil response. Modelling clearly revealed that the shallow substructure plays an important role in strong ground motion prediction using 3D simulation. The results show significant improvements in effective shaking duration and the peak ground velocity (PGV) distribution in terms of the accuracy achieved by our developed model. We recommend a high-resolution shallow substructure as an essential component in future seismic hazard analyses.


INTRODUCTION
Taiwan is located in a region of high seismicity on the boundary between the Philippine Sea Plate and Eurasian Plate.Of the ten disastrous Taiwan earthquakes in the past one hundred years (1898 -1997), five of these events occurred in the southwestern region (the 6 November 1904 Douliu earthquake, M L 6.1; 17 March 1906 Meishan earthquake, M L 7.1; 17 December 17 Zhongpu earthquake, M L 7.1; 5 December 1946 Xinhua earthquake, M L 6.1; and the 18 January 1964 Baihe Earthquake, M L 6.3).The southwestern plain is the largest alluvial plain on the island, comprising Tainan City and Chiayi County.Most researchers consider that there is a high probability that the next big earthquake in Taiwan will occur in this area.One analysis (Wen et al. 2014) has shown high earthquake event probabilities in the western foothills, which appear to be a highly active tectonic zone.As this region is faced with such a great earthquake threat, it is necessary to understand the strong earthquake motion characteristics there more clearly to reduce the potential seismic hazards.Many studies (Wen and Peng 1998;Lee et al. 2008a, b;Miksat et al. 2010) have indicated that the top soft soil layer can amplify ground shaking significantly during an earthquake.Furumura and Hayakawa (2007) conducted seismic-wave simulations using a 3D sedimentary structural model of Central Japan, and pointed out that long period ground motions with 7 s in Tokyo developed propagating through the sediments (3 -4 km) that overlie the bedrock.Takemura et al. (2015) conducted 3D finite-difference (FD) method simulations of the longperiod ground motions with 3D sedimentary velocity structure model for the Kanto Basin in Japan.They demonstrated that the lateral variation in about 1 km sediments determines the long-period (4 -8 s) ground motion characteristics in the northern Kanto Basin. Lee et al. (2008a, b) performed wave-field simulations within the Taipei basin using a (3D) discontinuous FD method.Their results indicated that the top soft sediments can amplify seismic waves and plays an important role when energy propagates in the Taipei basin.The inclusion of unconsolidated sedimentary layers in the subsurface structure is essential to improve the accuracy of numerical modelling and seismic hazard assessments.
Recently, many seismic tomographic studies (Wu et al. 2007(Wu et al. , 2009; Kuo-Chen et al. 2012;Huang et al. 2014) have provided higher-resolution 3D velocity images that incorporate an increasing number of stations and records in the Taiwan region.These large-scale velocity models make it possible to include heterogeneous velocity models in numerical modelling to improve the accuracy of wave propagation modelling (Miksat et al. 2010;Lee et al. 2014a, b).
In this study we attempt to construct a 3D velocity model of the western Taiwan plain, including the shallow substructure, with a low velocity and crustal structure, to improve 3D ground motion prediction for seismic hazard evaluation and risk management.However, there was no shallow substructure data available for inclusion in the 3D velocity model because the existing shallow velocity information is insufficient for regions other than the Taipei Basin.
The microtremor array technique is an effective tool that has been widely used to obtain the shear wave velocity structures in several areas in Taiwan (Satoh et al. 2001;Huang and Wu 2006;Kuo et al. 2009;Lin et al. 2009;Wu and Huang 2013).However, the non-uniform distribution of survey sites in the literature makes mapping the shallow velocity structures to these depths laborious and unsatisfactory.We recently processed array measurement data that were collected at forty-five uniformly distributed sites (Fig. 1) in the western plain of Taiwan between 2004 to 2007 and 2015.We used these measurements to estimate the S-wave velocity structure at each site (Kuo et al. 2016).We will then use these velocity structures as well as the six S-wave velocity structures obtained through large array measurements (Lin et al. 2009), which have a higher resolution to delineate the depth distribution of the engineering bedrock.
We construct 3D velocity models of the western plain of Taiwan with the inclusion of a shallow substructure.We also test and assess these models using numerical wave propagation modelling.Velocity field snapshots, synthetic waveforms, peak ground velocity amplification (PGV), and strong shaking duration are analysed.Our results point out the importance of the shallow substructure in numerical simulation applications.

SHALLOW SUBSTRUCTURE MODEL
The western coastal plain of Taiwan is low and flat.It is the largest plain on the island.This plain is filled with alluvium deposits, resulting in flat-lying Quaternary layers.In order to map the detailed shallow sub-layer structure, microtremor array measurements at forty-five sites were conducted in western Taiwan.The maximum likelihood frequency-wavenumber method was adopted to obtain the Rayleigh wave phase velocities as a function of frequency.A genetic algorithm technique based on the fundamental Rayleigh wave mode was then implemented to calculate the optimum S-wave velocity model at each site.The initial model results were then inverted using the conventional surface wave inversion method (Herrmann 1987) to remove dummy layers.The final result for each site was the best-fit model with the lowest variance.The velocity structure estimated using microtremor arrays mostly fit very well with the logging profiles (Engineering Geological Database for TSMIP, EGDT; Kuo et al. 2011Kuo et al. , 2012) ) at the same sites.From the Earthquake Engineering viewpoint the "Engineering Bedrock" definition is based on the following, to use the shallower underlying formation interfaced with S-wave velocity from 300 -700 m s -1 .The "Engineering Bedrock" definition is widely used in the strong motion prediction field.We assumed the S-wave velocity of the engineering bedrock to be 600 m s -1 and then used a gridding algorithm 'surface' in the Generic Mapping Tools software (Wessel and Smith 1998) for spatial interpolation to obtain the upper surface depths of the engineering bedrock.The Engineering Bedrock depth between 24 -383 m gradually increases from the eastern foot of the mountains to the western coastline.The depths of the Pliocene (Pan 1967(Pan , 1968) and upper Miocene (Shaw 1996) formation tops in the deeper part were remapped as our shallow sub-layer model and the average S-wave velocities were estimated (Lin et al. 2009).These velocities were then used to construct a 3D velocity model with a shallow substructure and a deep crustal model (Kuo-Chen et al. 2012;Huang et al. 2014) for seismic wave simulation and ground motion prediction using numerical modelling.The three main depth contours of the sub-layer were determined using Engineering Bedrock.The tops of the Pliocene and upper Miocene formations are displayed in Fig. 2. We created two model pairs (Fig. 3) for each crustal model with and without the shallow substructure in order to test and compare the simulation results.Model A is the original crustal model by Huang et al. (2014), the S-wave velocity of the cubic is from 0.65 -4.26 km s -1 , P-wave velocity is from 1.9 -7.4 km s -1 .This model considered borehole logging data to improve the shallow velocity structure.However, the velocity correction treated only the very near surface part due to the investigated borehole logging depth just approached about 30 m.We created Model A' by combining our shallow substructure model with Model A. The western part of the model filled three sub-layers with low velocity.The eastern part is a mountainous area that was not modified.The S-wave velocity of the cubic is from 0.4 -4.26 km s -1 , P-wave velocity is from 1.6 -7.4 km s -1 .Model B is the original crustal model by Kuo-Chen et al. (2012), the S-wave velocity of the cubic is from 2.5 -3.8 km s -1 , Pwave velocity is from 3.5 -7.7 km s -1 .We created Model B' by combining our shallow substructure model with Model B. The S-wave velocity of the cubic is from 0.4 -3.8 km s -1 , P-wave velocity is from 1.6 -7.7 km s -1 .Model B has higher velocity than Model A, especially near the surface.Even so, Model A has a little higher S-wave velocity than Model B in the deep part.Thus, Model A appears to have more velocity variation than Model B. We estimated the P-wave velocity of the engineering bedrock from PS-logging database (EGDT).The original P-wave velocity of Layers 2 and 3 are from seismic refraction surveys (Sato et al. 1970) and average by Lin et al. (2009).We gave slightly slower Pwave velocity 3.1 km s -1 for Layer 3 and separate two V P 2.1 and 1.8 km s -1 for Layer 2 according to our S-wave velocity profiles.The V P /V S ratios were calculated for each Layer.The V P /V S ratios are 4, 3, 1.75, and 1.72 from the surface to deep Layer.Thus, we think our estimated P-wave velocity values are reliable.So far the accurate 3D P-wave velocity structure is unavailable for this area.The velocities and Q values used in the model sub-layers are shown in Table 1.
The size of the model region is 100.9 × 109.9 km horizontally and 26.9 km vertically and it is located between 23 -24°N and 120 -121°E.A density model is also required for FD method simulation: we calculated it from P-wave velocities after Glaznev et al. (1996).The convert equation is given as below: when vp less than 5500 m s -1 : when vp large or equal 5500 m s -1 : ( ) ( ) .ln vp vp 1656 1068 3 18 The velocity and density values of the 3D model were interpolated using Matlab triangle-based cubic interpolation routine 'scatteredInterpolant' to become the grid points of the FD grid with dimensions 120 × 60 m in the horizontal and vertical directions, respectively.

NUMBERICAL MODELING
A 3D FD code (Furumura and Chen 2005; Furumura and Kennett 2005) for wave-field simulation was applied to simulate wave propagation in the western plain.Synthetic seismic wave propagation is calculated by solving the wave equations for motion using the 16 th -order staggered-grid FD method in the horizontal direction and a 4 th -order staggeredgrid FD method in the vertical direction, and 2 th -order in time.At the model sides and bottom damping (Cerjan et al. 1985) and a one-way absorbing boundary (Clayton and Engquist 1977) were applied to avoid reflection wave interference.The maximum frequency allowed by the 3D FD (Furumura and Chen 2005;Furumura and Kennett 2005) simulation was at least 1.1 Hz, according to the minimum S-wave velocity 0.4 km s -1 and the 0.12 km model grid size.The interval of each time step is 0.003 s.The modelling parameters are listed in Table 2.The total grid consisted of about 344 million points in one 3D model.We divided the model into 64 horizontal slices and conducted the simulation using parallel computing with communication by Message Passing Interface (MPI).The simulations were performed on the Institute of Earth Sciences, Academia Sinica supercomputer.The average computing time for a 60 s simulation is about 27 h.Parallel computation enables the large-scale computing required to efficiently assess and test these models.
We simulated the earthquake, M w = 6.1, that occurred beneath Chiayi City on 22 October 1999 at 23.53°E, 120.45°N with a focal depth of 16 km.The source focal mechanism is defined as strike, 356.6°; dip, 68°; and rake, 71.2° from CMT catalogue in the Broadband Array in Taiwan for Seismology (BATS) (http://bats.earth.sinica.edu.tw/), the preferred slip rupture model given by Chi and Dreger (2004).The event was located in the central region of our model.It was strong enough to trigger most stations in the dense strong-motion TSMIP network in this region (Taiwan Strong Motion Instrumentation Program; Tsai and Lee 2005).We simulated the event to compare observations and modelling results.We simulated wave propagation in each of the four different models (Fig. 3).Except for the velocity parameters, all modelling parameters input into the four different models were identical.In order to assess these models we compared the observed ground motions from this event with those simulated by our developed model and tomography model in terms of wave-fields, waveforms, PGV, and duration.The quality factor Q is an important parameter in the modelling of seismic-wave propagation to describe wave attenuation.The attenuation structure in sediments affected the amplitude of the surface wave packet (Takemura et al. 2015).We proposed a low velocity 3D structure as the sediment deposits for calculating the wave propagation in the western plain.Therefore, a suitable Q model assumption is required.We adopted the empirical relation between V S and Q values of Brocher (2008) to obtain the Q P and Q S models: We used an approximate Q model because the accurate values for this region were unavailable.

SIMULATION RESULT
We simulated the velocity waveforms for 86 stations that recorded the 22 October 1999 Chiayi earthquake.We selected ten stations of which eight were near the fault plane area and two stations were located on mountainous areas to validate our source assumption by comparing the synthetic and observed velocity waveforms.The distribution of 86 and 10 stations are shown in Fig. 4.Here we apply quantitative misfit criteria based on the correlation coefficient to evaluate the quality of our simulations: Where f is the observed waveform, g is the simulated waveform, T is the time length.Thus, when the observed and simulated waveforms are same, the misfit is 0. Conversely the misfit is 2. Figure 5 displays the simulated and observed waveforms obtained from Model A'.The velocity waveforms are band-pass filtered between 0.05 -0.3 Hz.Based on the fitted P-wave first arrival time, we shift the synthetic S-wave time in order to fit the observation.Most of the simulated and observed waveforms around the source area had good phase fit and had low misfit.Only few components had large misfit between the simulated and observed waveforms because the accuracy of our velocity model is not high enough.We think the source parameter is reliable in our FD method simulation although the simulated waves are still slightly faster than the observed.Developing a more accurate, higher resolution velocity model is essential.We did not shift the synthetic S-wave time except to validate our source assumption.All 3D synthetic velocity waveforms from Model A' and the observed velocity waveforms for sites corresponding to the TSMIP in the western plain of Taiwan are shown in Fig. 6.The solid squares denote the station locations.The synthetic and observed velocity waveforms are drawn together in red and blue, respectively, for comparison.The velocity waveforms were band-pass filtered from 0.02 -0.8 Hz. Figure 7 displays the waveform comparison between Model A' and the observed clearly for the sixteen stations labelled in Fig. 6.The waveforms are very complex and show long ground motion duration, especially in the southern part of the plain (e.g., CHY115, CHY044).Stations near the epicentre (e.g., CHY009, CHY046) show large amplitudes for body waves that are affected mainly by the source effect.Most of the synthetic waveforms reproduce the observed waveforms acceptably.

Vp (km s -1 ) Vs (km s -
Figure 8   reflections and mode conversions between the shallow sublayer and the free surface.Surface waves with large amplitudes appear particularly in the vertical component between 42 -60 s in the southern area.This feature can also be observed in real seismic waves.Conversely, Model A and B have a much shorter duration compared to Model A' and B', indicating that no surface wave is generated in these models.It should be mentioned that Model A exhibits longer shaking than Model B because of the low, very-thin velocity layer on its surface.Lee et al. (2015) proposed a new definition for the strong shaking duration called "effective shaking duration" (ESD).They first considered the time series window during which the amplitude is larger than 0.01 g, and then defined the ESD as the length of the time window for the accumulated energy within 5 -95% for the total energy.The total energy definition was given by McCann and Shah (1979).We adopt the definition by Lee et al. (2015) to calculate the ESD in order to quantify the fit of our developed models.The SAC (SEISMIC ANALYSIS CODE; Goldstein et al. 2003) differentiation routine was applied to obtain the acceleration waveforms to calculate the simulated ESD. Figure 9 shows a comparison between the simulated and observed ESD values.We calculated the PGV quadratic mean from all three components (band-pass filtered from 0.02 -0.8 Hz) as the peak value of the respective simulated or observed velocity waveforms.Figure 10 shows a comparison between the synthetic and observed PGV.The synthetic PGV maps are obtained from Model A' (Fig. 10a) and Model B' (Fig. 10b).Both PGV distributions show slightly overestimated amplitudes.Model A' has better PGV fit with observation than Model B'.However, they display similar patterns to their observed values.This seems to indicate that the PGV distributions are controlled mainly by the source radiation pattern, because the epicentre is located near the surface.Figure 11 shows PGV value and residual comparisons for Model A with A' and for Model B with B'.In order to evaluate the model's ability to predict ground motion, we did not apply any filters to the observed and synthetic waves in this section.For Model A and A' the lnErr v values of the horizontal components are similar at around 0.8.Any improvement in ground motion prediction from using the updated model is not apparent in this case.However, the lnErr v of the horizontal components reduces from 1.219 -0.883 for Model B and B'.PGV prediction has slightly improved because of the updated velocity model that includes shallow sub-layers in this case.The lnErr v of the vertical components are all around 0.5 -0.6 for the two model pairs.No improvement in PGV prediction in the vertical component due to the updated model is apparent.The ground motion prediction of the vertical component is merely adequate.Quality factor model or fault plane parameter simplification for the source may be responsible for its weak performance.The large misfit between the vertical component waveforms may also be a general numerical modelling problem (Lee et al. 2008b;Miksat et al. 2010).Nevertheless, the aim of this paper was accomplished.

DISCUSSIONS AND CONCLUSIONS
We simulated a M w = 6.1 earthquake that occurred near Chiayi City in our developed model as a test case to simulate 3D wave propagation to assess the influence of shallow substructures on ground motion.Although the synthetic waveforms obtained were still insufficiently accurate, the simulation was nonetheless valuable because it permitted assessing the influence of shallow substructures on ground motion.This in turn helped us to understand the seismic response associated with shallow sub-layers.We constructed 3D models that combine the shallow sedimentary substructure (Kuo et al. 2016) and the crustal and mantle structure of Taiwan (Kuo-Chen et al. 2012;Huang et al. 2014).The Topography effect of a local earthquake is negligible in sedimentary areas.Figure 10 in Lee et al. (2008a) shows the peak ground acceleration (PGA) amplification factor only increases at mountain tops and ridges for a local earthquake.The topography effect impacts the surrounding sedimentary area (Taipei basin) in the case of an earthquake occurs near the I-Lan region (Lee et al. 2009).In our study we did not include surface topography in the applied FD scheme.The source is located within sedimentary areas and we focused mainly on discussing the response of the western plain, which is characterized by flat unconsolidated deposits.Ignoring the topographic effect should be reasonable in this case.
The ESD was calculated in this study for the simulated accelerations using different models and observed records.The results show significant improvements in ESD in terms of the accuracy achieved by our developed model.Our model would be useful for 3D simulations for assessing seismic hazards, including building performance, potential landslides and so on.By comparing the simulation results from different models with corresponding observations, we calculated the PGV amplification factor due to the shallow  in Model A and A' has varied distribution.The northern areas have the largest amplification factors, reaching 300%, whereas the central area is around 50% or less than 50% (Fig. 12a).The mountainous areas have reduced ground motions which the largest amplification factors reaching -90% for both models.The largest amplification reduced area for Model A, A' is around northern part; for Model B, B' is near southern area.The phenomenon of reduced PGV amplification may due to more seismic energy being trapped in the low velocity substructure, causing less seismic energy propagating to the mountainous areas.Although the shallow substructure model is the same for both Models A' and B', the PGV amplification factors of these models are different.This indicates that the deeper velocity structure (> 5 km) still has an influence on PGV amplification.The velocity contrast between the source and surface is dominant in determining the degree of wave amplification and the level of velocity contrast between the soft layer and the basement affects surface waves generated by multiple reflections and mode conversions.It is worth noting that seismic wave amplification can be influenced by several factors, such as source radiation, topography, and substructure geometry.Each geometric factor could be frequency dependent.
We conclude that the shallow sublayers place high constraints on strong ground motions, especially by influencing the duration of strong shaking.Furthermore, the deep crustal model also influences the resultant ground motions.Therefore, investigations into high-resolution velocity structures and improved resolution of the deeper structure are necessary steps that will provide important information for the construction of more-accurate 3D models.

Fig. 1 .
Fig. 1.Distribution of the forty-five sites (yellow diamonds) for microtremor array measurements in the western coastal plain of Taiwan.Blue circles denote sites in the large microtremor arrays conducted by Lin et al. (2009).(Color online only)

Fig. 2 .
Fig. 2. Three major discontinuities in the western plain.(a) Depth contour of engineering bedrock (Layer1).(b) Depth contour of the top of the Pliocene formation (Layer2).(c) Depth contour of the top of the Miocene formation (Layer3).

Fig. 4 .
Fig. 4. Map of the Western plain with fault plane solution of BATS catalogue.The black rectangle is the fault dimension of the preferred slip model given by Chi and Dreger (2004).Diamonds are the TSMIP station distribution that recorded the 22 October 1999 Chiayi earthquake.Red diamonds are the stations we used to validate the source parameters.(Color online only) Fig.6.Site-by-site comparison between synthetic and observed velocity waveforms of the 22 October 1999 M w 6.1 Chiayi earthquake.The velocity waveforms obtained by Model A' are band-pass filtered between 0.02 -0.8 Hz.

Fig. 7 .
Fig. 7. Comparison between synthetic and observed velocity waveforms for sixteen stations.(Color online only) quantify the difference between the simulated and observed ESD.The observed and simulated ESD values for Model A and A' and their respective lnErr v values are displayed in Figs.9a and c. lnErr v is reduced from 0.632 -0.488.The observed and simulated ESD values for Model B and B' and their respective lnErr v values are displayed in Figs.9b and d. lnErr v is reduced from 0.755 -0.391.The results indicate that the inclusion of the shallow sub-layer in the models significantly improves the simulated ESD accuracy.

Fig. 9 .Fig. 10 .
Fig. 9. Comparison between observed ESD and synthetic ESD.(a) Correlation between ESD and hypo-distance for observed and synthetic ESD for Model A and A'; (b) correlation between ESD and hypo-distance for observed and synthetic ESDs for Model B and B'; (c) comparison between ESD residual and hypo-distance for Model A and A'; (d) comparison between ESD residual and hypo-distance for Model B and B'.(Color online only)

Fig. 11 .
Fig. 11.Comparison between observed and synthetic PGV for Model A, A', B, and B'.In this figure, the PGV values are not filtered.(Color online only)

Table 1 .
The velocity and Q value used in the model sub-layers.