3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses

TECTO-127230; No of Pages 18 Tectonophysics xxx (2016) xxx–xxx

Contents lists available at ScienceDirect

Tectonophysics journal homepage: www.elsevier.com/locate/tecto

3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses Dong Hun Lee a, Jung Mo Lee a,⁎, Hyun-Moo Cho b, Tae-Seob Kang c a b c

Department of Geology, Kyungpook National University, Daegu 41566, Republic of Korea Earthquake Research Center, Korea Institute of Geoscience and Mineral Resources, Daejeon 34132, Republic of Korea Department of Environmental Sciences, Pukyong National University, Busan 48513, Republic of Korea

a r t i c l e

i n f o

Article history: Received 31 July 2015 Received in revised form 21 August 2016 Accepted 30 August 2016 Available online xxxx Keywords: 3D velocity structure Receiver function Moho Gyeongju area of Korea

a b s t r a c t A temporary seismic array was in operation between October 2010 and March 2013 in the Gyeongju area of Korea. Teleseismic records of the seismic array appropriate for receiver function analysis were collected, and selected seismograms were split into five groups based on epicenters—the Banda-Molucca, Sumatra, Iran, Aleutian, and Vanuatu groups. 1D velocity structures beneath each seismic station were estimated by inverting the stacked receiver functions for possible groups. The inversion was done by applying a genetic algorithm, whereas surface wave dispersion data were used as constraints to avoid non-uniqueness in the inversion. The composite velocity structure was constructed by averaging the velocity structures weighted by the number of receiver functions used in stacking. The uncertainty analysis for the velocity structures showed that the average of 95% confidence intervals was ±0.1 km/s. The 3D velocity structure was modeled through interpolation of 1D composite velocity structures. Moho depths were determined in each composite velocity structure based on the AK135-F S-wave velocity model, and the depths were similar to the H-κ analysis results. The deepest Moho depth in the study area was found to be 31.9 km, and the shallowest, was 25.9 km. The Moho discontinuity dips in a southwestward direction beneath the area. A low velocity layer was also detected between 4 and 14 km depth. Adakitic intrusions and/or a high geothermal gradient appear to be the causes of this low velocity layer. The 3D velocity structure can be used to reliably assess seismic hazards in this area. © 2016 Elsevier B.V. All rights reserved.

1. Introduction Crustal velocity structure is an important research topic in seismology because various quantitative descriptions of earthquake-related phenomena include seismic velocity as a parameter. Velocity structures are essential in the determination of accurate hypocentral parameters from observed earthquake records. They are also useful for estimating ground motions due to earthquakes and explosions. Furthermore, accurate high-resolution data on crustal velocity structures improve the accuracy of earthquake reporting and the reliability of earthquake early warning systems. It allows seismologists to accurately determine focal parameters for setting up seismotectonics by correlating earthquakes with faults. An understanding of velocity structures also aids seismologists in reliably assessing seismic hazards. There are two major methods for exploring crustal velocity structures, namely crustal scale seismic profiling and earthquake data inversion. The former provides detailed (high frequency) absolute velocity structures using controlled sources, although costs are high ⁎ Corresponding author. E-mail addresses: [email protected] (D.H. Lee), [email protected] (J.M. Lee), [email protected] (H.-M. Cho), [email protected] (T.-S. Kang).

(e.g., Brown et al., 1987; Clowes et al., 1992; Kummerow et al., 2004). The latter provides deep regional (low frequency) velocity structures at low expense, although having an earthquake observation network is a prerequisite and results are limited by uncertainty due to tradeoffs between velocity and dimensions (e.g., Langston, 1979; Zhao et al., 1992; Lin et al., 2008). These methods complement each other. However, in Korea, there are difficulties in studying crustal velocity structures using local earthquake data because of low seismicity with relatively small magnitude earthquakes. In order to overcome these difficulties, the receiver function analysis technique was selected. The technique could utilize globally distributed large earthquakes and the study area was covered by a dense array of broadband seismic stations with the capability of recording earthquakes. The Korean Crust Research Team (KCRT) has performed crustal scale seismic experiments using controlled sources to investigate the 3D crustal velocity structures of the Korean Peninsula area (e.g., Cho et al., 2006; Kim et al., 2007; Cho et al., 2013). The KCRT2002 survey line was extended to the Gyeongju area (Cho et al., 2006); however, explosion of the planned source in this area was finally cancelled. The velocity cross-section of this area is therefore not as reliable as in other areas on the survey line. On the other hand, 1D velocity structures have been derived using receiver function analysis (Chang et al., 2004; Chang and

http://dx.doi.org/10.1016/j.tecto.2016.08.022 0040-1951/© 2016 Elsevier B.V. All rights reserved.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

2

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

Baag 2005, 2007; Yoo et al., 2007; Park et al., 2009; Han et al., 2010). Park et al. (2009) suggested that the Moho dip to the southwest beneath the HDB station in the Gyeongju area. Kang and Park (2006) derived tomograms for the southern Korean Peninsula area, showing velocity distributions of depth slice sections, even if these were of low resolution. Gyeongju has been the capital of the long-lasting Silla Dynasty, which has left many histories (57 BCE–935 CE). The historical record entitled “Samguksagi (History of three kingdoms)” reported information on abundant seismic activity in the Gyeongju area, and the event in March 779 CE which caused N 100 deaths, (Korea Meteorological Administration, 2012) appears to have been the one with the biggest impact. An increase in the number of large man-made structures due to industrialization and the discovery of Quaternary faults (e.g., Lee et al., 1999; Kee et al., 2007) in this area demand reliable seismic hazard assessment. A brief tectonic map of the Korean peninsula region is presented in the index map of Fig. 1 in order to help understand the regional tectonic setting. The purpose of this work is to estimate the 1D velocity structure beneath each station of the temporary seismic array in the Gyeongju area of Korea, and to model the 3D velocity structure through resampling and point kriging of the 1D velocity models. The final result is interpreted within the context of local geology and geophysical phenomena. 2. Receiver function and genetic algorithm Receiver function analysis is a common, well-established seismological method to estimate crustal velocity (e.g. Langston, 1979; Owens et al., 1984; Ammon et al., 1990; Chang and Baag, 2005; Yoo et al., 2007; Park et al., 2009). Teleseismic P waveforms contain information

on the source, the structure near the source and the receiver, and propagation effects (Cassidy, 1992). We applied the deconvolution technique proposed by Langston (1979) to process the data. This allows suppression of the earthquake source, the structure near the source, and the propagation effects. The receiver function therefore shows the relative response of the Earth structure under the receiver. The specific method is outlined in Langston (1979) and Ammon (1991). A genetic algorithm (GA) was used in the inversion; this is a stochastic global search method capable of near-optimal solutions for a multivariable function (Goldberg, 1989). The GA performs calculations with many random initial models within a prescribed range of the model parameters. The result does not strongly depend on the initial model. We chose GA as the inversion method because the result is weakly dependent on the initial model and we have little a priori information on crustal velocities. The uniform crossover method suggested by Syswerda (1989), in which each gene has the same probability, was adopted in this study. Receiver function analysis mainly utilizes traveltime differences between different phases and the velocity contrasts between layers; therefore, the inversion result is not unique (Ammon et al., 1990). To mitigate this tradeoff, surface wave dispersion data were used as additional, independent constraints. It is well known that receiver functions are sensitive to the shear wave velocity contrasts between layers; whereas surface wave dispersion is sensitive to the average shear wave velocities of layers. Given the high degree of complementarity of the information provided by these two types of data, their combination into collaborative inversion approaches is rather common (e.g., Ammon et al., 1990; Julia et al., 2000; Chang et al., 2004; Lawrence and Wiens, 2004; Park et al., 2009; Li et al., 2010; Shen et al., 2013). The phase velocities of Rayleigh waves obtained with the velocity model chosen by GA are compared with the observed phase velocities. If the difference is larger than 0.1-km/s, even at only one period, then the

Fig. 1. Location map of the temporary broadband seismic array of 18 stations, together with the permanent seismic station (HDB), in the Gyeongju area of Korea. The tectonic map shown in the inset is redrawn from Chough and Sohn (2010). Shaded areas in the inset map indicate Cretaceous sedimentary basins. Abbreviations are: IB; Imjingang Belt, SKTL; South Korean Tectonic Line, MTL; Median Tectonic Line.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx Table 1 Station information: number, code, coordinates, and brief description. Station Station Coordinates # Code Lat. (°N) Lon. (°E)

Elevation (m)

1

ONGT

35.5828 129.1628 96

2

GILC

35.5785 129.0936 135

3

EHWA

35.6613 129.3434 81

4

HTMJ

35.6545 129.2412 279

5

MIHO

35.6773 129.1719 159

6

SOHO

35.6700 129.0862 402

7

MALB

35.7487 129.3307 145

8

DUKC

35.7470 129.1784 78

9

NEIL

35.7493 129.0855 462

10

GAMP

35.8369 129.4957 84

11

YDON

35.8316 129.4132 134

12

DDON

35.8355 129.3236 185

13 14

HCCH BUSS

35.8064 129.1416 80 35.8396 129.0735 172

15

GALP

35.9071 129.4027 120

16

DJAW

35.8860 129.2898 92

17

BJTM

35.9198 129.2175 102

18

GSSW

35.8904 129.1609 94

19

HDB

35.7337 129.3991 146

Description

Surface bedrock Shale, Mudstone 14-m-deep drill well Alluvial deposit 8-m-deep drill well Granite at 4 m Surface bedrock Andesite Surface bedrock Andesite 14-m-deep drill well Granite at 7.5 m 22-m-deep drill well Soft rock 6-m-deep drill well Soft rock, Cobble stone at 4 m Surface bedrock Andesite Surface bedrock Weathered granite 22-m-deep drill well Tuff at 7 m 22-m-deep drill well Shale at 16 m Alluvial deposit Surface bedrock Tuff Surface bedrock Weathered breccia Surface bedrock Weathered conglom. and hornb. Surface bedrock Quartz porphyry Surface bedrock Shale, Mudstone Permanent station 65-m-deep drill well

Group 4

3

Table 2 List of earthquakes used in receiver function analyses. No.

Date (yyyy-mm-dd)

Time (UTC) (hh:mm:ss)

Lat. (°N)

Lon. (°E)

BAZ (°)

Dist. (°)

MW

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42

2010-10-08 2010-10-08 2010-10-08 2011-04-06 2011-06-13 2011-06-21 2011-06-24 2011-07-31 2011-08-16 2011-08-20 2011-08-20 2011-09-01 2011-09-02 2011-09-05 2011-10-23 2011-11-14 2011-12-13 2012-01-10 2012-02-02 2012-02-03 2012-02-05 2012-02-05 2012-03-09 2012-04-11 2012-04-11 2012-04-14 2012-04-15 2012-06-23 2012-07-06 2012-07-25 2012-08-10 2012-08-11 2012-08-18 2012-08-26 2012-09-05 2012-09-26 2012-10-08 2012-10-17 2012-10-20 2012-12-02 2012-12-11 2012-12-21

05:43:08 03:49:10 03:26:13 14:01:43 14:31:22 02:04:15 03:09:39 14:34:47 11:03:57 16:55:02 18:19:23 06:14:38 10:55:53 17:55:11 10:41:22 04:05:11 07:52:11 18:36:59 13:34:40 03:46:21 00:15:38 16:40:39 07:09:50 10:43:10 08:38:36 22:05:26 05:57:40 04:34:53 02:28:22 00:27:45 18:37:43 12:23:18 09:41:52 15:05:37 13:09:10 23:39:55 11:43:31 04:42:30 23:00:32 00:54:22 06:18:27 22:28:08

2.831 51.287 51.374 1.615 2.515 −11.479 52.050 −17.016 −2.331 −18.365 −18.311 −12.357 52.171 2.965 38.722 −0.949 0.041 2.433 −17.827 −17.378 −18.894 −17.948 −19.125 0.802 2.327 −18.972 2.581 3.009 −14.657 2.707 52.633 38.329 −1.315 2.190 −12.476 51.592 −4.472 4.232 −13.552 −16.975 0.533 −14.344

128.217 −175.180 −175.361 97.097 126.457 165.551 −171.836 171.579 128.002 168.143 168.218 166.657 −171.708 97.893 43.513 126.910 123.030 93.210 167.133 167.277 168.917 167.226 169.613 92.463 93.063 168.741 90.269 97.896 167.340 96.045 −167.421 46.826 120.096 126.837 166.513 −178.295 129.129 124.520 166.564 167.645 126.231 167.286

181.88 50.55 50.42 228.54 185.07 136.90 49.35 135.16 182.00 138.94 138.84 136.44 49.17 228.81 301.21 183.90 190.53 233.43 139.49 139.09 138.60 139.48 138.15 232.77 233.49 138.80 236.51 228.85 137.32 230.65 48.28 299.54 194.99 184.34 136.66 50.04 180.16 188.98 137.35 138.51 185.22 137.17

33.17 42.57 42.46 45.72 33.57 58.86 44.65 66.83 38.36 66.01 66.00 60.21 44.73 44.15 66.48 37.02 36.57 47.57 65.01 64.72 66.87 65.17 67.44 49.30 47.75 66.83 49.44 44.12 62.49 45.52 47.34 64.31 38.31 33.87 60.22 40.61 40.49 32.05 61.26 64.58 35.58 62.20

6.2 6.0 6.4 6.0 6.3 6.0 7.3 6.1 6.1 7.2 7.1 6.0 6.9 6.7 7.1 6.3 6.0 7.2 7.1 6.1 6.1 6.1 6.7 8.2 8.6 6.2 6.2 6.1 6.3 6.4 6.2 6.4 6.3 6.6 6.0 6.4 6.1 6.0 6.2 6.1 6.0 6.7

velocity model is discarded. This method was used to resolve the nonuniqueness problem of the receiver function analysis. The effectiveness of receiver function analysis using GA and constrained by surface wave dispersion has been numerically shown by Chang et al. (2004), and the

Group 3 3.8

Group 5 Group 2 Group 1

Velocity (km/s)

3.7

3.6

3.5

3.4

3.3 0 Fig. 2. Locations of teleseismic events used in this work. The events are categorized into five groups – Banda-Molucca (Group 1), Sumatra (Group 2), Iran (Group 3), Aleutian (Group 4), and Vanuatu (Group 5) – based on epicentral distributions. Lines represent the great circle paths linking each event and the study area.

10

20

30

Period (second) Fig. 3. Rayleigh wave phase velocities at the HDB station used as constraints in receiver function inversions (redrawn from Chang and Baag, 2005).

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

4

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

(a)

(b)

(c)

(d)

Fig. 4. Estimated velocity structures beneath stations (a) BJTM, (b) BUSS, (c) DDON, and (d) DJAW. The stacked and synthesized receiver functions (blue line and red line, respectively) for each group are on the right side of the velocity structure. The blue arrow indicates the low velocity layer, while the red arrow indicates the identified Moho discontinuity. Thin red dotted lines present the limits of the prescribed velocity range. Estimated velocity structures beneath (e) DUKC, (f) EHWA, (g) GALP, and (h) GAMP. Estimated velocity structures beneath (i) GILC, (j) GSSW, (k) HCCH, and (l) HTMJ. Estimated velocity structures beneath (m) MALB, (n) MIHO, (o) NEIL, and (p) ONGT. Estimated velocity structures beneath (q) SOHO, (r) YDON, and (s) HDB.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

5

(e)

(f)

(g)

(h)

Fig. 4 (continued).

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

6

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

(i)

(j)

(k)

(l)

Fig. 4 (continued).

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

7

(m)

(n)

(o)

(p)

Fig. 4 (continued).

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

8

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

(q)

(r)

(s)

Fig. 4 (continued).

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

technique has been successfully applied in numerous works (e.g., Chang and Baag, 2005; Park et al., 2009; Han et al., 2010).

9

using data from the HDB station by Chang and Baag (2005). Fig. 3 illustrates each phase velocity. Assuming the Poisson solid, the density of pffiffiffi each layer was assigned by ρi ¼ 0:32 3 βi þ 0:77 (Berteussen, 1977).

3. Data processing 4. Results Teleseismic events recorded by a seismic array consisting of 18 stations in the Gyeongju area of Korea were used in this study. The array was installed in a roughly 10-km-spaced grid pattern, consisting of five points in the east-west direction and four in the north-south direction, and in operation between October 2010 and March 2013. Table 1 gives station numbers, codes, coordinates, and a brief description of each station, while Fig. 1 shows station locations. Ten seismometers were installed on surface bedrock, seven in 8-inch drill wells, and one on surface alluvium. The HDB downhole permanent station in this area was included in the analysis. Table 1 also gives details of drill well depths and local geology. We selected 42 events with magnitudes larger than 6 and with epicentral distances between 30° and 90°. A 10-s window was chosen before and after P wave arrival time to compute signal-to-noise ratio (SNR). The average SNR was 7.5 dB. The events were categorized into five groups, with events within a single group located within a range of 10° in back azimuth and epicentral distance. The five categories were the Banda-Molucca group, the Sumatra group, the Iran group, the Aleutian group, and the Vanuatu group, as shown in Fig. 2. Table 2 provides the list of the earthquakes used in receiver function analyses. The radial and transverse components were computed through rotation, using the back azimuths for ground surface seismometers. The data recorded by eight seismometers installed in the drill wells (including HDB) were rotated to the angle at which the horizontal P phase matches the vertical P phase at the teleseismic frequency band. This rotation angle neglected horizontal inhomogeneity, and transverse receiver functions sensitive to the dipping velocity discontinuity were not analyzed in this work. The water-level parameter was set to 0.01 and the Gaussian filter parameter was set to 2.5, corresponding to an effective 0.5 Hz cutoff frequency. Receiver functions were stacked to enhance SNR (Cassidy, 1992). Our grouping criteria of 10° for back azimuth and epicentral distance were tighter than the stacking bounds suggested by Owens et al. (1984). The quality of individual receiver functions were examined and divided into subgroups based on waveform similarity criteria, and the major subgroup with high quality and similarity was chosen for stacking. We used Rayleigh wave phase velocities for the period between 10 and 25 s, with a 5-s interval estimated

4.1. 1D velocity structure beneath each station The 1D velocity structure beneath each seismic station was estimated by inversion of the stacked receiver functions constrained by surface wave dispersion for the five groups, and the composite velocity structure was constructed by averaging the velocity structures weighted by the number of receiver functions in stacking. The Moho depth was determined based on the AK135-F model (Kennett et al., 1995; Montagner and Kennett, 1996), which displays an S-wave velocity of 3.9 km/s at the base of the crust. The Moho discontinuity was assumed to correspond to the point in the composite velocity structure where the S-wave velocity jump was largest within the depth range with an Swave velocity higher than 3.9 km/s (Fig. 4). Instead, the depth of the second largest velocity jump closest to the largest one was chosen as the Moho depth if the difference between the two velocity jumps was b0.05 km/s and if the depth of the second largest jump was preferable in comparison with Moho depths at nearby stations. Station elevations were not corrected for comparison purposes, and station elevations for the 1D velocity structures were thus fixed to 0 km. The composite velocity structures were then resampled at 0.1 km depth intervals and the station elevations were corrected. The corrected velocity structures were used for comparison beneath stations and for modeling the 3D velocity structure. The inverted models at most stations displayed low velocity layers (blue arrows in Fig. 4) with variable thicknesses between 4 and 14 km depth in the upper crust, while the Moho was located at depths between 26 km and 32 km (red arrows in Fig. 4). A summary of results from the receiver function analysis was given in Table 3, including the number of receiver functions, the depth range of the low velocity layers, and the elevation corrected Moho depth beneath each station. The uncertainty of the composite velocity models beneath 19 stations was analyzed utilizing the jackknife resampling scheme (Tukey, 1958; Efron and Stein, 1981). 95% confidence intervals were evaluated and presented with the composite velocity models in Fig. 5, and the numerical values were listed in Table S1 in the supplementary material.

Table 3 Summary of results from receiver function analyses. Station number

Station code

Number of individual RFs

Low velocity layer depth range (km)

Moho depth from Vs; DVS (km)

Moho depth from H-κ analysis; DH ­κ (km)

Moho depth difference; (km) DVS − DH ­κ

VP/VS ratio

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

BJTM BUSS DDON DJAW DUKC EHWA GALP GAMP GILC GSSW HCCH HTMJ MALB MIHO NEIL ONGT SOHO YDON HDB

34 26 13 15 15 28 22 20 10 14 10 28 26 20 23 23 13 19 16

5–8 5–12 6–10 5–12 5–8 6–8 5–12 6–14 5–12 6–8 6–12 4–8 6–10 5–12 6–10 4–8 6–10 6–12 5–8

27.90 31.80 27.85 27.90 27.95 27.90 27.90 25.90 31.90 31.90 29.90 29.70 27.90 27.85 29.55 31.90 31.60 27.90 27.90

28.50±2.1 30.50±1.4 27.50±3.3 24.40±1.1 27.45±2.7 25.95±1.4 26.90±2.0 24.40±3.0 30.30±2.8 31.90±3.4 29.90±3.0 27.20±1.1 26.90±1.5 29.35±1.8 28.55±1.4 31.40±2.8 28.10±2.2 24.40±0.9 24.90±2.4

−0.60 1.30 0.35 3.50 0.50 1.95 1.00 1.50 1.60 0.00 0.00 2.50 1.00 −1.50 1.00 0.50 3.50 3.50 3.00

1.76±0.08 1.81±0.02 1.78±0.20 1.71±0.02 1.83±0.03 1.81±0.02 1.80±0.03 1.84±0.07 1.79±0.03 1.66±0.09 1.78±0.03 1.80±0.02 1.78±0.02 1.74±0.02 1.86±0.02 1.75±0.03 1.86±0.06 1.89±0.01 1.85±0.05

Moho depths are station elevation corrected values.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

10

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

0

0

0

0

10

10

10

10

20

20

20

20

30

30

30

30

BJTM 40

2

BUSS 3

4

5

40

2

4

5

40

2

3

GAMP 4

5

40

0

0

0

0

10

10

10

10

20

20

20

20

30

30

30

30

GILC 40

Depth (km)

EHWA 3

2

HTMJ 3

4

5

40

2

MALB 3

4

5

40

2

4

5

40

0

0

0

10

10

10

10

20

20

20

20

30

30

30

30

NEIL 2

ONGT 3

4

5

40

2

3

5

40

2

3

5

40

0

0

0

10

10

10

10

20

20

20

20

30

30

30

30

DDON 2

3

DJAW 4

5

40

2

4

5

40

0

0

0

10

10

10

20

20

20

30

30

30

GSSW 40

2

3

2

HCCH 4

5

40

2

2

DUKC 3

3

3

5

40

2

5

3

4

5

3

4

5

3

4

5

GALP 4

5

3

40

2

95% C.I. Comp. Vel.

SOHO 4

4

HDB 4

0

40

2

YDON 4

3

MIHO 3

0

40

2

4

5

S-wave velocity (km/s) Fig. 5. The 95% confidence intervals of the composite velocity models for 19 stations estimated using the jackknife resampling scheme. The lower and upper limits are shown in blue lines, and the composite velocity is the thick red line.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

The composite velocity uncertainty was generally about ±0.1 km/s and did not exceed ±0.2 km/s. The characteristic features of the velocity structure beneath each station that were not presented in Table 3 were visually summarized in Fig. 6. The velocity structures fell into two categories, one that shows smoothly varying velocity layers and one that fluctuated. The former

11

category could further be divided into two groups containing gradually increasing velocity layers in the lower crust. An intracrustal velocity discontinuity was also sometimes apparent in the former category. To visualize the spatial variation of the 1D velocity models, the composite velocity structure beneath each station was plotted on the location map of broadband seismic stations (Fig. 7). As shown in Fig. 7, the

Fig. 6. Characteristic features of velocity structures beneath each station not already presented in Table 3. The black horizontal line indicates the intracrustal velocity discontinuity. Abbreviations are: LVL; low velocity layers, DLVL; distinctive low velocity layers, SVVL; smoothly varying velocity layers, FVL; fluctuating velocity layers, and GIVL; gradually increasing velocity layers.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

12

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

Fig. 7. 1D composite velocity structures beneath each station with the location map of seismic stations. Blue and red arrows indicate the depths of the low velocity layer and the Moho discontinuity identified from the S-wave velocities, respectively.

Moho depth decreased eastwards. We also applied the H-κ analysis proposed by Zhu and Kanamori (2000) in order to estimate the error bounds of the obtained Moho depth values. The results were presented in Fig. 8 and detailed values of Moho depths with standard deviations and VP/VS ratios were listed in Table 3. The Moho depth determined by the S-wave velocity was compared to the one determined by H-κ analysis for each station and also presented in Table 3. In general, the trends showed good agreement. The differences of six stations (DJAW,

EHWA, HTMJ, SOHO, YDON, and HDB) exceed the standard deviation obtained by H-κ analysis. For EHWA and HTMJ, the Moho depths were determined as 28 and 30 km, respectively, by choosing the largest velocity jump, and were 2 km deeper than those determined by H-κ analysis. If the second largest jump was chosen, the depths were 26 and 28 km, respectively, and agreed well with those obtained through H-κ analysis. Based on the criteria for the Moho depth by S-wave velocity, the depths of the second largest jump at the HDB station and the largest jump at the

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

13

Fig. 8. Amplitude distributions in the H-κ domain of the grid search at each station.

SOHO station were chosen. If those of the largest jump of HDB and the second largest jump of SOHO were chosen, they agreed well with those from H-κ analysis. To attempt to explain the disagreement between results from the two analysis methods, we reviewed standard deviations from H-κ analysis, but found no obvious explanation. The transverse components of receiver functions, which were sensitive to horizontal inhomogeneity, were also reviewed. Large transverse

components were found in receiver functions from five stations (BJTM, EHWA, HCCH, MIHO, and ONGT). 4.2. 3D velocity structure Depths of the composite velocity models of the 19 stations were corrected using station elevations with a resolution of 50 m. The

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

14

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

elevation-corrected Moho depths were interpolated using 1-km interval point kriging in the N-S and E-W directions. The resulting Moho depth distribution was presented in Fig. 9(a). The Moho depth distribution using 14 station results (excluding BJTM, EHWA, HCCH, MIHO, and ONGT stations, which had large transverse components in the receiver functions) was obtained for comparison purposes, and was presented in Fig. 9(b). The Moho depth distribution using H-κ analysis was also presented in Fig. 9(c). Standard deviations of H-κ analysis results were presented in Fig. 9(d). The general trends of the three obtained Moho depth distributions were similar to one another, excluding small deviations. Zhu and Kanamori (2000) confirmed that receiver functions with H-κ analysis could provide “a very good point measurement of the crustal thickness under a broadband station”. The standard deviation of the H-κ analysis was considered to be a good measure of model parameter uncertainty. Any unknown factor that could reduce reliability, such as the VP/VS ratio, was avoided in the interpretation of receiver function analysis results, and the confidence interval of the velocity model was narrow (±0.1 km/s in average). The Moho depth distribution determined from S-wave velocities was expected to be similar to the one obtained with H-κ analysis.

Elevation-corrected composite velocities for every 100 m slice were interpolated using 1-km interval point kriging in the N-S and E-W directions. We combined the 2D velocity slices for each depth and produced the 3D velocity structure. Fig. 10 presented the velocity distributions of depth slice sections at 3 km depth intervals, enabling visualization of 3D velocity structures. To enhance visualization, the E-W and N-S vertical slice sections for every 6-km interval were also shown in Figs. 11 and 12, respectively, where the Moho discontinuity was shown by dotted lines. In general, the Moho discontinuity gradually deepened from east to west. As shown in Fig. 11(e) and (f), and Fig. 12(b), (c), (d), and (e), the low velocity layer was evident around the 6–10 km depth interval in the north of the study area. 5. Discussion and conclusions A seismic array in the Gyeongju area of Korea, comprising 19 stations, was in operation from October 2010 to March 2013 (with the exception of the permanent HDB station). The receiver function analysis technique was applied to model the 3D velocity structure of this area. Teleseismic records appropriate for receiver function analysis were

Fig. 9. Moho depth and standard deviation distributions of the study area: (a) Moho depth distribution using the results of all 19 stations from S-wave velocities, (b) Moho depth distribution from S-wave velocities excluding the results of 5 stations with large transverse components, (c) Moho depth distribution from H-κ analysis for all 19 stations, and (d) standard deviation distribution of depths from H-κ analysis.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

15

2 km

5 km

8 km

11 km

14 km

17 km

20 km

23 km

26 km

29 km

32 km

35 km

Fig. 10. S-wave velocity distributions for (a) 2, (b) 5, (c) 8, (d) 11, (e) 14, (f) 17, (g) 20, (h) 23, (i) 26, (j) 29, (k) 32, and (l) 35 km depth slice sections. The contour interval is 0.02 km/s, with a velocity color scale. Seismic station locations are shown in (a) to enable precise positioning.

collected from these stations and selected seismograms were split into five groups, based on epicenters the Banda-Molucca, Sumatra, Iran, Aleutian, and Vanuatu groups. For each station, receiver functions of all available records of all possible groups were computed. The 1D velocity structures beneath each seismic station were estimated by inverting the stacked receiver functions using a genetic algorithm,

and surface wave dispersion was used as a constraint. The composite velocity structure beneath each station was constructed by averaging the velocity structures, weighted by the number of receiver functions in stacking since SNR increased with stacking number and the inversion result generally became more reliable. The velocity model parameter uncertainty was analyzed using the jackknife

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

16

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

Fig. 11. Velocity distributions of E-W vertical slice sections at (a) 3941, (b) 3947, (c) 3953, (d) 3959, (e) 3965, and (f) 3971 N-S UTM km. The contour interval is 0.1 km/s, with a velocity color scale. The Moho discontinuity is shown by dotted lines.

Fig. 12. Velocity distributions on N-S vertical slice sections at (a) 511, (b) 517, (c) 523, (d) 529, (e) 535, and (f) 541 E-W UTM km. The contour interval is 0.1 km/s, with a velocity color scale. The Moho discontinuity is shown by dotted lines.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

17

Fig. 13. Geothermal gradient map from the New and Renewable Energy Data Center (2012). The seismic station locations are shown to enable precise positioning.

scheme, and the confidence interval was ± 0.1 km/s on average. The 3D velocity structure was modeled through point kriging of 1D composite velocity structures, and the Moho depth beneath each station was derived (Table 3). The first major result of this work concerned the Moho depth distribution. This generally increased in a southwest direction, although there were local variations. A southwestward dipping Moho was reasonable, considering that the Korean Peninsula was located close to the eastern margin of the Eurasian plate. The deepest Moho depth was found to be 32 km while the shallowest was 26 km. The reliability of the Moho depths determined by S-wave velocities was examined by comparing the results with those of H-κ analysis. The results also appeared to be reasonable by comparison with previous independent results; the local depth variation of 6 km within a square with 40 km sides seemed to be large. The Moho depth beneath the HDB station was determined as 28 km, which was compatible with previous estimates, including the southwestward dipping 28-km deep Moho (Park et al., 2009) and the westward dipping 28-km deep Moho (Chang and Baag, 2005) determined for the same station. A second major result of this work concerned the low velocity layer that was found over most of the study area at a depth of 4–14 km. 3D seismic tomograms of Kang and Park (2006) showed low velocity layers in this area; however, the horizontal cell sizes (40 km × 40 km) of their tomograms are too large to be compared in detail. Two alternative explanations are proposed for the presence of the low velocity layers. First, it is hypothesized that the low velocity layer reflects mantle upwelling in the slab window of the subduction zone (Gerya et al., 2006; Qi et al., 2007). Late-Cretaceous I-type granite such as hornblende biotite granodiorite has been reported in the eastern part of the study area by Ko et al. (1996). These authors analyzed the major, trace, and rare earth element compositions of the granitic rocks. Ishihara and Jin (2005) suggested a modified definition of adakite based on the composition of granite satisfying the following conditions: SiO2 (N 56%), Al2O3 (N 15%), Na2O (N3.5%), and Sr/Y (N20). The I-type granite in the study area yields SiO2, Al2O3, and Na2O ratios of 65.7%, 15.17%, and 2.63%, respectively. The major and trace elements in this granite were mapped into the adakite region by Kiji et al. (2000). The plot of Sr/Y vs Y ratios, known as the discrimination diagram for the adakite, also gave convincing evidence of adakitic rocks. The I-type granite of the study area classified by Ko et al. (1996) therefore corresponds to adakitic rocks. This high-alumina, volcanic rock demonstrates a link between subduction

of young oceanic crust and production of intermediate to felsic igneous rocks from partial melting (Thorkelson and Katrin Breitsprecher, 2005). Late-Cretaceous adakitic rocks, reported in the southwestern of Japan, are usually formed due to the melting of young and hot oceanic crust at a subduction zone (Takahasi et al., 2005). An alternative explanation for the low velocity layer is the high geothermal gradient. The area where the low velocity layer is distinct (north and west portions of the study area) corresponds to the area where the geothermal gradient is high (Fig. 13), based on mapping conducted by the New and Renewable Energy Data Center (2012). High temperatures cause thermal expansion, reducing seismic velocity due to lower density. The results of this study may also be the result of a combination of these two factors. The 3D velocity structure modeled in this study can be applied to the improvement of seismic hazard assessment reliability in the Gyeongju area. Supplementary data to this article can be found online at http://dx. doi.org/10.1016/j.tecto.2016.08.022. Acknowledgements This work was funded by the Korea Meteorological Administration Research and Development Program under Grant KMIPA 20157060. The data used in this work were collected by the “Establishment of probabilistic assessment basis for seismic safety of NPP site” project of Korea Hydro and Nuclear Power Corporation, Ltd. The authors are grateful to Dr. V. Sallares, Managing Editor of the Deep Seismix S.I., two anonymous reviewers, and Prof. C.-E. Baag for their help in improving the manuscript, and Prof. S.-J. Chang for his help with the research. References Ammon, C.J., 1991. The isolation of receiver effects from teleseismic P waveforms. Bull. Seism. Soc. Am. 81, 2504–2510. Ammon, C.J., Randall, G., Zandt, G., 1990. On the nonuniqueness of receiver function inversions. J. Geophys. Res. 95, 15303–15318. Berteussen, K.A., 1977. Moho depth determinations based on spectral ratio analysis of NORSAR long-period P waves. Phys. Earth Planet. Inter. 15, 13–27. Brown, L., Wille, D., Zheng, L., DeVoogd, B., Mayer, J., Hearn, T., Sanford, W., Caruso, C., Zhu, T.-F., Nelson, D., Potter, C., Hauser, E., Klemperer, S., Kaufman, S., Oliver, J., 1987. COCORP: new perspectives on the deep crust. Geophys. J. R. Astr. Soc. 89, 47–54. Cassidy, J.F., 1992. Numerical experiments in broadband receiver function analysis. Bull. Seism. Soc. Am. 82, 1453–1474.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022

18

D.H. Lee et al. / Tectonophysics xxx (2016) xxx–xxx

Chang, S.-J., Baag, C.-E., 2005. Crustal structure in Southern Korea from joint analysis of teleseismic receiver functions and surface-wave dispersion. Bull. Seism. Soc. Am. 95, 1516–1534. Chang, S.-J., Baag, C.-E., 2007. Moho depth and crustal VP/VS variation in southern Korea from teleseismic receiver functions: implication for tectonic affinity between the Korean Peninsula and China. Bull. Seism. Soc. Am. 97, 1621–1631. Chang, S.-J., Baag, C.-E., Langston, C.A., 2004. Joint analysis of teleseismic receiver functions and surface wave dispersion using the genetic algorithm. Bull. Seism. Soc. Am. 94, 691–704. Cho, H.M., Baag, C.-E., Lee, J.M., Moon, W.M., Jung, H., Kim, K.Y., 2006. Crustal velocity structure across the southern Korean Peninsula from seismic refraction survey. Geophy. Res. Lett. 33, L06307. http://dx.doi.org/10.1029/2005GL025145. Cho, H.M., Baag, C.-E., Lee, J.M., Moon, W.M., Jung, H., Kim, K.Y., 2013. P- and S- wave velocity model along crustal scale refraction and wide-angle reflection profile in the southern Korean peninsula. Tectonophysics 582, 84–100. Chough, S.K., Sohn, Y.K., 2010. Tectonic and sedimentary evolution of a cretaceous continental arc-backarc system in the Korean peninsula: new view. Earth Sci. Rev. 101, 225–249. Clowes, R.M., Cook, F.A., Green, A.G., Keen, C.E., Ludden, J.N., Percival, J.A., Quinlan, G.M., West, G.F., 1992. Lithoprobe: new perspectives on crustal evolution. Can. J. Earth Sci. 29, 1813–1864. Efron, B., Stein, C., 1981. The jackknife estimate of variance. Ann. Stat. 9, 586–596. Gerya, T.V., Connolly, J.A.D., Yuen, D.A., Gorczyk, W., Capel, A.M., 2006. Seismic implications of mantle wedge plumes. Phys. Earth Planet. Inter. 156, 59–74. Goldberg, D.E., 1989. Genetic Algorithms in Search, Optimization and Machine Learning. Addison-Wesley. Han, S.I., Lee, J.M., Kang, T.S., 2010. 1D crustal velocity structure beneath broadband seismic stations in the Okcheon Fold Belt of Korea by receiver function analysis. Geoscience J. 14, 57–66. Ishihara, S., Jin, M.S., 2005. The Jecheon granite-a western extension of the Hida granite? Proceedings of Spring Meeting. Kor. Soc. of Econ. Environ. Geol. 16-17. Julia, J., Ammon, C.J., Herrmann, R.B., Correig, A.M., 2000. Joint inversion of receiver function and surface wave dispersion observations. Geophys. J. Int. 143, 99–112. Kang, I.B., Park, J.H., 2006. Application of seismic tomography to the region in and near southern Korean Peninsula. Econ. Environ. Geol. 39, 507–524 (in Korean with English abstract). Kee, W.S., Kim, B.C., Hwang, J.H., Song, K.Y., Kihm, Y.H., 2007. Structural characteristics of Quaternary reverse faulting on the Eupcheon Fault, SE Korea. J. Geol. Soc. Korea 43, 311–333. Kennett, B.L.N., Engdahl, E.R., Buland, R., 1995. Constraints on seismic velocities in the earth from travel times. Geophys. J. Int. 122, 108–124. Kiji, M., Ozawa, H., Murata, M., 2000. Cretaceous adakitic Tamba granitoids in northern Kyoto, San'yo belt, southwest Japan. Jpn. Mag. Mineral. Petrol. Sci. 29, 136–149. Kim, K.Y., Lee, J.M., Moon, W., Baag, C.E., Jung, H., Hong, M.H., 2007. Crustal structure of the Korean Peninsula from seismic waves generated by large explosions in 2002 and 2004. Pure Appl. Geophys. 164, 97–113. Ko, J.S., Yoon, S.H., Lee, S.W., 1996. Petrology and geochemical characteristics of A-type granite with particular reference to the Namsan Granite, Kyeongju. J. Petrol. Soc. Korea 5, 142–160 (in Korean with English abstract). Korea Meteorological Administration, 2012. Historical Earthquake Records in Korea (2–1904). KMA press, Seoul in Korean.

Kummerow, J., Kind, R., Oncken, O., Giese, P., Ryberg, T., Wylegalla, K., Scherbaum, F., TRANSALP Working Group, 2004. A natural and controlled source seismic profile through the Eastern Alps: TRANSALP. Earth Planet. Sci. Lett. 225, 115–129. Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J. Geophys. Res. 84, 4749–4762. Lawrence, J.F., Wiens, D.A., 2004. Combined receiver-function and surface wave phase-velocity inversion using a niching genetic algorithm: application to Patagonia. Bull. Seism. Soc. Am. 94, 977–987. Lee, B.J., Ryoo, C.-R., Chwae, U., 1999. Quaternary faults in the Yangnam area, Kyongju, Korea. J. Geol. Soc. Korea 35, 1–14 (in Korean with English abstract). Li, H., Bernardi, F., Michelini, A., 2010. Surface wave dispersion measurements from ambient seismic noise analysis in Italy. Geophys. J. Int. 180, 1242–1252. Lin, F., Moschetti, M.P., Ritzwoller, M.H., 2008. Surface wave tomography of the western United States from ambient seismic noise: Rayleigh and Love wave phase velocity maps. Geophys. J. Int. http://dx.doi.org/10.1111/j.1365-246X.2008.03720.x. Montagner, J.P., Kennett, B.L.N., 1996. How to reconcile bodywave and normal-mode reference Earth models? Geophys. J. Int. 125, 229–248. New and Renewable Energy Data Center, 2012. Geothermal Gradient Map. http://kredc.kier.re. kr/kier/. Owens, T.J., Zandt, G., Taylor, S.R., 1984. Seismic evidence for an ancient rift beneath the Cumberland Plateau, Tennessee: a detailed analysis of broadband teleseismic P waveforms. J. Geophys. Res. 89, 7783–7795. Park, S.J., Lee, J.M., Ryu, I.-C., 2009. 1D velocity structure beneath broadband seismic stations in the Cretaceous Gyeongsang Basin of Korea by receiver function analyses. Tectonophysics 472, 158–168. Qi, C., Zhao, D., Chen, Y., 2007. Search for deep slab segments under Alaska. Phys. Earth Planet. Int. 165, 68–82. Shen, W., Ritzwoller, M.H., Schulte-Pelkum, V., Lin, F.-C., 2013. Joint inversion of surface wave dispersion and receiver functions: a Bayesian Monte-Carlo approach. Geophys. J. Int. 192, 807–836. Syswerda, G., 1989. Uniform crossover in genetic algorithms. In: Schaffer, J. (Ed.), Proceedings of the Third International Conference on Genetic Algorithms. Morgan Kaufmann Publishers, Los Altos, CA, pp. 2–9. Takahasi, Y., Kagashima, S.-I., Mikoshiba, M., 2005. Geochemistry of adakitic quartz diorite in the Yamizo Mountains, central Japan: implications for Early Cretaceous adakitic magmatism in the inner zone of southwest Japan. Island Arc 14, 150–164. Thorkelson, D.J., Katrin Breitsprecher, K., 2005. Partial melting of slab window margins: genesis of adakitic and non-adakitic magmas. Lithos 79, 25–41. Tukey, J.W., 1958. Bias and confidence in not quite large samples. Ann. Math. Stat. 29, 614–623. Yoo, H.J., Herrmann, R.B., Cho, K.H., Lee, K., 2007. Imaging the three-dimensional crust of the Korean peninsula by joint inversion of surface wave dispersion and teleseismic receiver functions. Bull. Seism. Soc. Am. 97, 1002–1011. Zhao, D., Hasegawa, A., Horiuchi, S., 1992. Tomographic imaging of P and S wave velocity structure beneath northeastern Japan. J. Geophys Res. 97, 19909–19928. Zhu, L., Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver functions. J. Geophys. Res. 105, 2969–2980.

Please cite this article as: Lee, D.H., et al., 3D crustal velocity structure beneath the broadband seismic array in the Gyeongju area of Korea by receiver function analyses, Tectonophysics (2016), http://dx.doi.org/10.1016/j.tecto.2016.08.022