- Email: [email protected]
Crustal structure in the Changbaishan volcanic area, China, determined by modeling receiver functions
Tectonophysics 386 (2004) 157 – 175 www.elsevier.com/locate/tecto
Crustal structure in the Changbaishan volcanic area, China, determined by modeling receiver functions E.A. Hetlanda,*, F.T. Wub, J.L Songb,1 a
Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139-4307, USA b Department of Geological Sciences and Environmental Studies, State University of New York at Binghamton, Binghamton, NY 13902, USA Received 21 July 2003; accepted 2 June 2004 Available online 22 July 2004
Abstract During 1998–1999, we installed a temporary broadband seismic network in the Changbaishan volcanic region, NE China. We estimated crustal structure using teleseismic seismograms collected at the network. We detected a near surface region of strong anisotropy directly under the main volcanic edifice of the volcanic area. We modeled 109 receiver functions from 19 broadband stations using three techniques. First we used a bslant-stackingQ method to model the principal crustal P reverberation phases to estimate crustal thickness and the average crustal P to S speed ratio (v p/v s), assuming an average P-wave velocity in the crust. We then estimated crustal S-wave velocity (v s) and v p/v s profiles by modeling stacked receiver functions using a direct search. Finally, we inverted several receiver functions recorded at stations closest to the main volcanic edifice using least squares to estimate v s velocity profiles, assuming a v p/v s value. The results from the three estimation techniques were consistent, and generally we found that the receiver functions constrained estimates of changes in wave speeds better than absolute values. We resolved that the crust is 30–39 km thick under the volcanic region and 28–32 km thick away from the volcanic region, with a midcrust velocity transition at about 10–15 km depth. We estimated that the average crust P-wave velocity is about 6.0–6.2 km/s surrounding the main volcanic region, while it is slightly lower in the vicinity of the main volcanic edifice. The estimates of v p/v s were more ambiguous, but we inferred that the bulk crustal Poisson’s ratio (which is related to v p/v s) ranges between 0.20 and 0.30, with a suggestion that the Poisson’s ratio is lower under the central volcanic region compared to the surrounding areas. We resolved low S-wave velocities (down to about 3 km/s) in the middle crust in the region of the main volcanic edifice. The low velocity anomaly extends from about 5–10 to 15–25 km below the surface, probably indicating a region of elevated temperatures. We were unable to determine if partial melt is present with the data we considered in this paper. D 2004 Elsevier B.V. All rights reserved. Keywords: Changbai Mountains; PASSCAL; Crustal structure; Receiver functions; Tianchi volcano
* Corresponding author. E-mail address: [email protected] (E.A. Hetland). 1 Now at US Geological Survey, 384 Woods Hole Road, Woods Hole, MA 02543 USA. 0040-1951/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2004.06.001
158
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
1. Introduction The Changbaishan volcanic region is located in Jilin Province in northeastern China, along the border between China and North Korea (Fig. 1). There are two main volcanic edifices in the Changbaishan region, the Wangtian and the Tianchi, which are of basaltic and andesitic to rhyolitic composition underlain by a large Cenozoic basaltic shield (Liu, 1983; Liu et al., 1998). The Wangtian was last active 2.12 Ma and has approximately 400 km2 of exposed volcanic rock associated with it (Fan et al., 1999). The Tianchi volcano became active about 1.66–1.48 Ma years ago (Jin and Wei, 1994) and is still active, with the most recent activity in 1903 (Cui et al., 1995). The Tianchi eruption of 1000AD was one of the largest recorded in China in the past 2000 years (Liu et al., 1997). Until recently, little geophysical research has been conducted in the Changbaishan region. Regional Swave tomography by Guo et al., (1996) indicated a low velocity anomaly at 25–75 km depth below the Tianchi volcano (also see Liu et al., 1998). Additionally, a magnetotelluric study done by Tang (1998) revealed a region of low resistivity below 10–15 km depth extending from the Tianchi volcanic crater to about 30 km to the north. The region of low resistivity might indicate a region of partial melt or a region of high fluid content, but was interpreted to be most probably a magma chamber by Tang (1998). The Moho depth was determined to be about 30–35 km by modeling surface and Moho reflected phases from a regional earthquake (Shin and Baag, 2000). Based on geochemical analysis, Shangguan and Sun (1997) have determined that the gases released from local hot springs are from the top of a residual, mantle-derived magma chamber in the crust, located possibly as shallow as 5 km below the Tianchi crater. The water from hot springs in the Changbaishan region are about 80.75 8C and have been determined to be meteoric water circulating through regions of hot rock (Shangguan et al., 1997). In the summer of 1998, the State University of New York at Binghamton, in collaboration with the Research Center of Geophysical Exploration (RCGE) of the Chinese Seismological Bureau, installed 19 portable, broadband seismic stations in NE China, in
the region of the Changbaishan volcanic area (Table 1 and Fig. 1). The stations were borrowed from IRIS PASSCAL and each station comprised a Guralp 3T sensor and a six-channel Ref Tek recorder. We established stations at 18 sites in late June 1998. Eight sites were occupied until late September 1998, while 10 sites were occupied until March 1999. The one remaining station was installed in September 1998 and removed in March 1999. In August 1998, the RCGE conducted an active source, wide-angle reflection/refraction experiment, and analysis of profiles near the Tianchi volcano indicated a low Pwave velocity anomaly under the Tianchi volcano at about 10–20 km depth, and resolved fairly uniform Moho depths of about 30–35 km away from the volcanic center, thickening to about 40 km under the Tianchi edifice (Zhang et al., 2002; Song et al., in preparation). In this study, we present a determination of crustal structure using teleseismic receiver functions throughout the Changbaishan region. The receiver function method has been applied to many tectonic problems (e.g., Owens and Zandt, 1997; Julia` et al., 1998; van der Hilst et al., 1998; Chevrot and van der Hilst, 2000; Julia´ et al., 2000; Zhu and Kanamori, 2000) and volcanic regions (e.g., Langston, 1979; Darbyshire et al., 2000; Du and Foulger, 2001).
2. Receiver functions 2.1. Source equalization The first 20 or so seconds following direct P in a teleseismic seismogram principally contains rupture effects, source-side structure, and receiver-side structure. The source-side signal results from near source reflected P waves and S to P conversions. Receiverside signal consists of P to S conversions and P and S reflections near the receiver, principally in the crust (Fig. 2). Direct P is incident on the Moho nearly vertically for a teleseismic distance source. Therefore, the amplitude of S conversions are larger on the horizontal components, while the P signals are larger on the vertical component. Since the horizontal components of the seismogram principally contain receiver-side effects, it is possible to equalize the source effects and retain the receiver-side effects by
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
159
Fig. 1. Map of the temporary broadband station locations, circle type indicates the approximate deployment of each station. The Tianchi and Wangtian volcanic craters are indicated by white triangles, rivers, lakes and ocean are black, and the border between China and North Korea is a thick dashed line.
deconvolving the vertical component from the radial component leaving a radial receiver function (e.g., Phinney, 1964; Langston, 1979; Ammon, 1991). In
order to avoid division by small numbers, the equalization procedure of Langston (1979), modified by Ammon (1991), uses a water-level deconvolution
160
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
Table 1 Locations of the temporary broadband stations, the geographic region of each station relative to the Tianchi volcano (bTianchiQ signifies the station is near the main volcanic edifice; see Fig. 13), and the number of receiver functions (No.) determined at each station used in this study Station
Latitude (8N)
Longitude (8E)
Region
No.
Station
Latitude (8N)
Longitude (8N)
Region
No.
BACH CBAI DAND DRAG JANA JIYU LING PSIA TILL YABA
42.514 42.058 40.146 42.504 41.362 42.379 41.799 42.950 42.267 43.003
128.158 128.063 124.391 129.031 126.512 126.795 126.957 126.078 123.914 129.489
Tianchi Tianchi west north south south south west west north
5 1 8 5 19 5 3 8 3 7
CANY CH2A DHAA FROG JILA LIAO MANG THAA WUSU
41.951 43.801 43.344 41.420 43.728 41.272 41.948 41.697 42.326
128.005 125.450 128.198 128.174 126.664 123.284 127.594 125.998 127.277
Tianchi west north south west west south south south
3 3 16 7 3 5 3 3 2
where signal with spectral power below some threshold, called the water level, is discarded. This equalization method includes a Gaussian filter in the deconvolution, which acts as a low-pass filter removing noise in the seismograms. Because the receiver function method uses crustal converted and reflected phases to determine crustal structure, the region of crust sampled depends on the offset of these crustal reverberations (Fig. 2). The offsets of the reverberations depend on the ray parameter ( p) and the crustal thickness (or more generally the depth to the deepest interface). We estimated p from the event depth and range using the IASP91 model of P global seismic travel times (Kennett and Engdahl, 1991). The lateral extent from the station that the crust is sampled by a receiver function is at least the distance to the penetration point of the first bounce P arrival (Fig. 2). In this region where the crust is expected to be about 30–35 km thick, with an average P wave velocity of 6.0–6.5 km/s, and for the data set we used, the receiver
functions sample the Moho about 40–50 km from the stations. Furthermore, Sheriff and Geldart (1982) have demonstrated that the lateral resolution of seismically imaged interfaces is approximated by the radius of the Fresnel zone, and as the radius of the Fresnel zone increases, the P wave samples a larger area of an interface. Since the radius of the Fresnel zone increases with depth, the lateral resolution over an interface decreases with depth. In this area, the lateral resolution of P on the Moho is about 10–15 km. Energy propagating off the back azimuth will decrease the amplitude of the SV waves on the radial component and will not be resolved on the equalized radial receiver function (e.g., Ammon, 1991; Zhang and Langston, 1995). For this reason it is crucial to rotate the seismogram accurately, and strong lateral heterogeneity and receiver-side anisotropy will complicate this. In practice, equalization is not perfect, thus we considered all the equalized receiver functions to be possibly corrupted with noise, and we only
Fig. 2. Cartoon showing direct P (labeled Pp) and the dominant crustal reverberations, where, except for mantle P, lowercase and uppercase indicate upgoing and downgoing rays, respectively (a), and the associate receiver function (b); after Ammon et al. (1990).
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
tried to model dominant arrivals and arrivals which were consistent for all traces.
161
band seismometers were located south of the network, since the highest concentration of seismicity meeting the above criterion was in the seismic belt from Indonesia through Papua New Guinea to the Fiji islands (Fig. 3). We used one event from the Mariana Trench. We found only three events from the north (the Aleutian Trench region) and one from the west (located in Iran). Due to equipment failures, power outages, periods of unusually high anthropogenic noise levels, and differing durations of deployment not all of the events were recorded at all of the stations. In this study, we used 109 seismograms recorded at the 19 stations (Table 1). We verified the orientations of all stations, as well as the back azimuths of all events, by rotating the unequalized horizontal seismograms and visually inspecting the P energy on the radial and tangential
2.2. Data In this study, out of 41 earthquakes considered, we used 29 events below 48 km depth (one event was shallow, with a depth listed as 33 km) and at a range of 30–908 (Table 2 and Fig. 3). The first criterion was used to avoid complex near-source structure effects, while the second was used to ensure that the direct P arrived well before secondary P phases. We also used only large events to ensure that the signal was strong. The minimum magnitude event we found to equalize well was 5.3 Mb, although this depended on the signal-to-noise ratio of a given seismogram. Most of the accepted events detected by the temporary broad-
Table 2 Earthquakes used in this study Year
Day
Time
Latitude (8N)
Longitude (8E)
Depth (km)
Magnitude
Location
1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1999 1999 1999 1999 1999 1999 1999 1999 1999 1999
190 190 197 206 231 240 245 245 248 251 258 263 264 271 284 300 319 329 361 019 028 028 035 036 037 039 061 063 075
14:45:39.9 19:39:43.9 11:56:36.4 02:39:23.3 17:41:43.6 12:40:58.7 08:37:29.9 18:52:42.2 17:09:31.4 09:10:03.0 08:35:51.5 21:21:55.1 06:52:41.1 13.34:30.4 12:04:54.7 21:16:21.0 02:44:12.3 18:05:25.7 00:38:26.7 03:35:33.8 08:10:05.4 18:24:25.2 19:28:00.8 11:39:45.1 21:47:59.4 16:22:54.3 21:15:30.0 05:38:26.5 03:32:16.2
30.487 60.530 11.040 13.608 6.179 0.154 5.410 29.695 4.051 13.257 5.624 7.770 0.262 8.194 21.040 2.921 21.589 7.859 21.632 4.596 52.886 4.582 4.033 12.616 12.853 7.908 51.870 28.343 2.683
178.994 153.219 166.160 166.867 147.745 125.018 126.764 178.790 152.041 144.007 151.637 106.951 122.467 112.413 179.110 128.623 176.504 158.622 176.376 153.235 169.123 153.658 95.283 166.966 166.697 147.734 179.430 57.193 125.674
130.0 145.0 110.0 44.0 55.0 66.0 50.0 230.0 192.0 141.0 83.0 67.0 147.0 152.0 624.0 61.0 149.0 48.0 144.0 114.0 67.0 101.0 56.0 213.0 90.0 70.0 72.2 33.0 114.0
6.4ML 6.2ML 6.9ML 6.2ML 5.3MB 6.2MB 6.8ML 5.4MB 5.5MB 5.8MB 6.2ML 5.5MB 6.1MB 6.5ML 5.9ML 5.9MB 6.2ML 6.2ML 6.6ML 6.5ML 6.6ML 6.3ML 5.9ML 5.9ML 7.3MS 5.9ML 5.5MB 6.5MS 5.6MB
KERMADEC ISLANDS, NEW ZEALAND SOUTHERN ALASKA SANTA CRUZ ISLANDS VAUATU ISLANDS EASTERN NEW GUINEA REG., P.N.G. SOUTHERN MOLUCCA SEA MINDANAO, PHILIPPINE ISLANDS KERMADEC ISLANDS, NEW ZEALAND NEW BRITAIN REGION, P.N.G. MARIANA ISLANDS NEW BRITAIN REGION, P.N.G. JAWA, INDONESIA MINAHASA PENINSULA, SULAWESI JAWA, INDONESIA FIJI ISLANDS REGION HALMAHERA, INDONESIA FIJI ISLANDS REGION SOLOMON ISLANDS FIJI ISLANDS REGION NEW IRELAND REGION, P.N.G. FOX ISLANDS, ALEUTIAN ISLANDS NRE IRELAND REGION, P.N.G. NORTHERN SUMATERIA, INDONESIA SANTA CRUZ ISLANDS SANTA CRUZ ISLANDS EASTERN NEW GUINEA REG., P.N.G. A RAT ISLANDS, ALEUTIAN ISLANDS SOUTHERN IRAN TALAUD ISLANDS, INDONESIA
Note that the Celebes Sea region includes the Philippine Islands and Molucca Sea events.
162
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
Fig. 3. Epicentral locations (circles) of the earthquakes used in this study and great circle path (dashed lines) from the epicenter to the seismic array (approximated as the location of the Tianchi volcano). The square indicates the area of Fig. 1, and the black lines are major plate boundaries.
components. For a typical seismogram, when rotated correctly, all of the energy of direct P is on the radial component until about 5 s later when the Ps phase (S conversion at the Moho) is detected (Fig. 4a). We were unable to rotate all of the energy off the tangential component of seismograms recorded at CBAI, where there is equal amplitude energy on both the tangential and radial components in the first several seconds (Fig. 4b). Moreover, this energy appears to be out of phase on the two components, which is typical of anisotropy observations. Only seismograms of events to the south of CBAI were affected, and this led us to conclude that there is a region of possible anisotropy in the crust to the south of CBAI. We observed only a small hint of this phenomenon on seismograms from CANY, located about 10 km to the south of CBAI, and since the seismic waves to both of these stations sampled
similar lower crust but more distinct upper crust, we surmised that the region of anisotropy is in the shallow crust under the Tianchi crater. The anisotropy may be the result of cracked rock corresponding to the geothermal activity here, and it may be coincident with the source of gases detected in hot springs in this area. Due to this complexity in the seismograms, we only used one seismogram recorded at CBAI, which was the record of an event, located to the north, that was not affected by the anisotropy. We do not present an analysis of the region of potential anisotropy in this paper. 2.3. Analysis After identifying all usable seismograms, we determined the receiver functions using a constant Gaussian width parameter of 2.4 Hz and varying
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
163
Fig. 4. Radial and tangential components of seismograms of two events located in the Santa Cruz islands recorded at two stations: (a) FROG, located about 70 km to the south of the Tianchi crater and (b) CBAI, located adjacent and slightly to the north of the Tianchi. Time is relative to the approximate arrival time of the initial P wave, and events are to the south of the stations.
water levels (from 510-5 to 10-2). We determined the water level largely by the quality of the original seismogram and stability of the equalization. We interpreted the radial receiver functions using three analysis techniques which estimate differing sets of model parameters: a bslant-stackingQ procedure, a direct search forward modeling scheme, and a leastsquares inversion (LS). In all of these methods, we assumed that the crust is locally homogeneous, isotropic, and that the Moho is not dipping. Since the receiver functions sample relatively localized regions of the crust, we only assumed the crust to be homogeneous over these scales. We could not account for anisotropy in these analyses, but we did not use any seismograms that showed signs of
anisotropy. We used the results of the three methods to cross-check each other. While using the above three techniques does not provide independent measurements, due to the differing methods and model parameters estimated, the three methods indicate robustness of the final results. 2.3.1. Slant-Stacking The slant-stacking method has the ability to resolve simultaneously crustal thickness (H) and the P to S wave speed ratio (v p/v s), thereby partially resolving the thickness/velocity tradeoff (Zhu and Kanamori, 2000; Chevrot and van der Hilst, 2000). In this calculation, only three travel times, relative to direct P, are needed: Ps (P–S conversion at the Moho, where,
164
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
except for mantle P, lowercase and uppercase indicate upgoing and downgoing rays, respectively, in the crust; Fig. 2), PpPs, and PpSs (or PsPs). Following the
arguments of Chevrot and van der Hilst (2000), we weighted each of these arrivals the same. Chevrot and van der Hilst (2000) did not use the PpSs phase, but
Fig. 5. (a) Receiver functions determined at CANY station, with the predicted travel time curves estimated from the optimum solution of H and v p/v s determined by slant-stack analysis (v¯ p=6.0 km/s). (b) Solution surface constructed by a direct search over H and v p/v s with the optimum solution (star) and the 0.9 contour line. Scale is normalized and all solutions with negative sum are white.
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
we found including it better constrained the solution with these data. The differential travel-time equations depend on H, v p/v s, the average crustal P-wave
165
velocity (v¯ p) and the ray parameter ( p). The ray parameter is determined from the event depth and range, and we assumed values of v¯ p of 6.0 and 6.5 km/s.
Fig. 6. (a) Receiver functions determined at JANA station, with the predicted travel time curves estimated from the optimum solution of H and v p/v s determined by slant-stack analysis (v¯ p=6.5 km/s). (b) Solution surface constructed by a direct search over H and v p/v s with the optimum solution (star) and the 0.9 contour line. Scale is normalized and all solutions with negative sum are white.
166
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
Poisson’s ratio (m) can be determined from v p/v s through the relation ! 1 1 1 m¼ 2 2 vp =vs 1 and is a useful quantity for crustal studies since m is sensitive to mineral composition (e.g., Christensen, 1996), the presence of cracks and fluid content (e.g., O’Connell and Budiansky, 1974), and partial melt (e.g., Mavko, 1980). At each station, for a given choice of H and v p/v s, we summed the amplitudes of the receiver functions, at the predicted differential travel times of the Ps, PpPs, and PpSs phases: X s H; vp =vs ¼ aðtPPs Þ þ a tPPpPs a tPPpSs R
where R is a set of receiver functions and a(t) is the amplitude of a receiver function at time t. We grid searched over (H,v p/v s) parameter space, and we took the optimum solution to be at the maximum of the solution surface (Figs. 5 and 6). In the case of a
sufficient number of receiver functions and clear arrivals, the model space will show one global maximum (Fig. 6). In the neighborhood of the maximum, the parameter surface (s) has a trend which is related to the v p/v s–H tradeoff. Previously Zhu and Kanamori (2000) used the curvature of s in the two parameter directions to infer errors on the estimated parameters. Noting that this method does not account for the anticorrelation between v p/v s and H, Chevrot and van der Hilst (2000) used a boot-strapping method, where they systematically removed several traces and determined the solution for each realization, then taking the standard deviation of these solutions to estimate the errors. While this method of error estimation accounts for the obvious anticorrelation, we found that in this study it underestimated the errors, possibly due to the relatively small number of traces we had for any given station, the data quality, or the equalization method we used. We determined that the boot-strapping method underestimated the errors through forward trials, where we perturbed the optimum solution until the change in travel-time curves was beyond the resolution of the data. We
Fig. 7. Estimated crustal thickness (H, top) and Poisson’s ratio (v, bottom) from the slant-stacking analysis, using v¯ p=6.0 and 6.5 km/s, as a function of distance from the Tianchi crater. Error bars are the minimum and maximum values of H or v of the error ellipses (e.g., Figs. 5 and 6).
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
found that the 90% contour line of s was an effective proxy for the errors (Figs. 5 and 6). Since information on v¯ p is scant in this area, we performed slant stacking using both 6.0 and 6.5 km/s. While the values of H and m changed according to the value of v¯ p we used, the trends of changes in H and m across the broadband network remained consistent and the values were consistent with neighboring stations (Fig. 7). Due to the small number of receiver functions determined at CBAI and CH2A, we did not include these stations in the slant-stacking analysis. 2.3.2. Direct search We also modeled the observed receiver functions with a direct search algorithm, which uses a nearest-
167
neighbor search method (Sambridge, 1999). The nearest-neighbor algorithm directs a Monte Carlo search to a global solution, and in that it is similar to a genetic algorithm or a simulated annealing algorithm. The nearest-neighbor algorithm (NA) has been shown to be effective for complex problems and has been applied successfully to receiver function modeling, solving for v p/v s, layer thicknesses and v s at the top and bottom of each model layer (Sambridge, 1999). As opposed to other direct search algorithms, the NA is only dependent on one tuning parameter. Due to complexity of search space and number of parameters per layer, in practice we were unable to get good results using more than four or at most five model layers. Thus, the NA tended to pick up the
Fig. 8. Stack of all of the receiver functions determined at CANY station from events in the Celebes Sea region (thin line, top panel) with the best fit model (dashed line). The best fit model (solid line, bottom panels) was found using the NA direct search. Fit is normalized such that the best fitting model has misfit of 1.0, and the all of the models with at least twice the misfit of the best model are shown (white to gray lines) along with the average (dashed line) of the shown models.
168
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
major velocity discontinuities, while resolving velocity gradients between the discontinuities (Fig. 8). Solving for v p/v s is a large asset to this method. We found the NA to be more sensitive to changes in v p/ v s, which indicate horizons of changing properties, and less sensitive to the actual value of v p/v s (e.g., Fig. 9). We inverted stacked receiver functions using the NA. By stacking several receiver functions with similar back azimuths and ray parameters, it was possible to increase the signal-to-noise ratio of receiver functions. Due to the small number of receiver functions determined at most stations, we could not get a stacked receiver function for LING, MANG, TILL, or WUSU. At TILL station we modeled a single receiver function which had a particularly high signal-to-noise ratio. At most of the other stations there were enough receiver func-
tions to get at least one stacked function to use in this analysis, which we modeled to determine velocity profiles (Fig. 8). We determined three velocity profiles from three stations on the outskirts of the volcanic region (BRIT, DHAA, and JANA; Fig. 9), and two velocity profiles for stations DAND, JIYU, and YABA. Except for the event used at TILL, all events we used in this phase of the study were from south to southeast directions. So that resulting velocity profiles did not contain a discontinuity when the data did not require one, we searched over a large range of layer thicknesses and specified the same range of velocities at the top and bottom of adjacent layers. While we were unable to model the entire receiver function using a model with 4–5 layers, we were able to model the dominant arrivals (Figs. 8 and 9). In all cases, the Moho depths determined with the NA direct search were consistent
Fig. 9. Stacked receiver functions modeled using a NA direct search. Receiver functions were determined at JANA station and are from the New Britain Islands, Celebes Sea and Fiji Islands regions. Caption is as in Fig. 8.
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
169
Fig. 10. Comparison of the crustal thickness (H) and average v p (v¯ p) estimated at CANY and JANA stations by the slant-stack and NA direct search methods. Diamond and black crosses are the H and v¯ p determined by the NA method (diamond represents the one estimate at CANY, and size of cross indicates standard deviation of the three estimates at JANA). Gray circles are H estimated in the slant-stack analysis for the indicated v¯ p, error bars are shown as in Fig. 7.
with those obtained in slant-stacking, if the differences in average crustal v p were accounted for (Fig. 10). Additionally, at the stations where we were able to determine more than one velocity model, the results were consistent (Fig. 9).
2.3.3. Least squares The last technique we used to model the receiver functions was the least-squares inversion of Ammon et al. (1990). The least-squares inversion (LS) method requires an initial velocity model composed of an
Fig. 11. Two receiver functions (thin lines, left panels) determined at CANY station from events in the Celebes Sea region on days 240 (cele240) and 245 (cele245) of 1998, and the best fit modeled receiver functions determined through LS inversion (dashed thick lines, right panels) using a starting model composed of 40 one kilometer layers (initial velocity profile is thin gray line). The velocity models determined from modeling the stack of cele240 and cele245 with the NA direct search (Fig. 8) are shown for reference.
170
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
arbitrary number of layers. The v s in each layer is then estimated, holding the layer thicknesses constant and assuming a constant Poisson’s ratio of 0.25 (Ammon et al., 1990). We needed to apply smoothing constraints in the inversion in order to avoid checker-boarding in
final solution. The LS technique is more sensitive to velocity discontinuities than to the absolute magnitude of the velocities (Ammon et al., 1990). Our motivation for using this inverse technique was to determine higher resolution crustal velocity
Fig. 12. Velocity models determined by LS inversion of the receiver functions at stations closest to the Tianchi volcanic crater. Individual velocity models are shown as thin gray lines, and the average of all velocity profiles for each station is shown with thick lines. Gray lines projecting from the stations show the approximate extent from the stations that each receiver function sampled the crust along the great circle path. Cross-hatched regions show the approximate regions of the crust sampled by the receiver functions which detected a midcrust low-velocity anomaly. The locations of the Tianchi and Wangtian volcanic craters (triangles) and the Chinese national border (dashed line) are shown for reference.
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
profiles than the NA direct search allowed. We used a starting model with 40–45 1-km-thick layers. We applied varying levels of smoothing until we were able to model the dominant phases while avoiding
171
overcomplexity and checker-boarding in the resultant model. Fig. 11 shows the results of inverting two receiver functions determined from events in the Celebes Sea region recorded at CANY station, where
Fig. 13. Velocity models determined from modeling receiver functions recorded by the Changbaishan broadband network, grouped according to region. Average velocity models, across the indicated stations, determined using the NA direct search (smooth bold lines). The average velocity models from each NA direct search (solid, thin gray lines) are shown to indicate the spread. The average velocity models from the 40–45 layer LS inversion (stepped lines) are shown for reference. NA velocity models from BACH and CANY are the results of single waveform analysis. Gray lines projecting from the stations show the approximate extent from the stations that each receiver function sampled the crust along the great circle path, and right and left leaning striped backgrounds indicate v p/v s as in Fig. 14a and b, respectively.
172
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
3. Discussion
volcanic region, and 30–39 km in the volcanic region, with the thickest crust closest to the Tianchi crater (Figs. 7 and 13). All of the events recorded at CANY station are from the south to southeast, so the region of Moho sampled is to the south and southeast. Thus, receiver functions recorded at CANY station constrain crustal thickness south of the Tianchi crater, under the Wangtian volcanic edifice. Using direct search and least squares inversion results, we resolved the thickest crust in the area extending from FROG north to BACH, including possibly MANG (Figs. 12 and 13). We determined the v¯ p in the regions surrounding the main volcanic area (the Tianchi and Wangtian craters) to be about 6.0–6.2 km/s. The lowest v¯ p we found using the stacked waveforms were for BACH, CANY and YABA stations (all below 6.0 km/s with YABA being the lowest at v¯ pc5.7 km/s). Furthermore, all of the velocity profiles determined for stacked receiver functions contained a midcrust boundary at about 10–15 km depth, with the exception of YABA station, whose velocity profiles indicated a midcrust transition at 20–25 km depth (Fig. 13).
3.1. Crustal thickness and bulk velocity structure
3.2. Midcrust low velocity anomaly
Using the slant-stacking algorithm, we found the Moho to be at about 28–32 km depth away from the
The velocity profiles we obtained from modeling stacked receiver functions from BACH and CANY
the receiver function from the event on day 245 has a higher signal-to-noise ratio than the event on day 240. We only performed these higher resolution inversions using traces recorded at stations closest to the center of the volcanic region (BACH, CANY, DRAG, FROG, MANG, WUSU and one trace recorded at CBAI from the north, which was not as affected by anisotropy as the records from the south were). The crustal velocity structure of this region was of most interest and we wished to constrain the results of the NA direct search. We inverted a total of 26 traces recorded at these seven stations. In these inversions, we needed to apply relatively large smoothing (20– 40%) because of the overall data quality and noise level. There were five traces for which we could not obtain stable or consistent results and we discarded them. Furthermore, we only obtained one velocity profile at CBAI and WUSU stations. The resulting velocity profiles were consistent for similar areas of crust and are consistent with those determined above using the NA method (Fig. 12).
Fig. 14. Profiles of v p/v s determined through modeling stacked receiver functions through NA direct search. (a) Average v p/v s profiles (gray lines) determined at DAND, DRAG, DHAA, FROG, JANA, THAA and YABA stations (corresponding to the right leaning striped stations in Fig. 13) and the average of the set (thick line). (b) Average v p/v s profiles determined at LIAO, JILA, JIYU and PSIA stations (left leaning striped) and the average of the set (thick line). (c) Average v p/v s profiles determined from modeling stacked receiver functions at BACH and CANY.
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
both indicate a pronounced high velocity layer in the upper crust, with a low velocity anomaly in the middle crust (Fig. 12). The high velocity layer starts at about 5–10 km depth, and is about 10 km thick at BACH and 10–15 km thick at CANY. The S-wave velocities in the upper layer are over 4 km/s, while the velocities decrease to about 3 km/s or slightly lower in the middle crust. The high velocity layer is more pronounced in the velocity profile at CANY than at BACH, and the stacked receiver functions from both stations sampled crust to the south of the stations. We further investigated the crust in the volcanic region by modeling single receiver functions with least-squares inversions (Ammon et al., 1990). The velocity profiles indicated an ubiquitous low velocity zone in the midcrust which appears to be localized in the region extending from FROG to BACH station, possibly further to the north (Fig. 12). The low velocity anomaly is thicker and deeper south of CANY than under the north slope of the Tianchi volcano. Receiver functions that sampled crust near FROG station resolved a slight low velocity region at greater depth in the crust (Fig. 12). Receiver functions recorded at MANG do not suggest a low velocity anomaly and we concluded that the western extent of the low velocity anomaly is located under the western flank of the Changbaishan (Fig. 12). One velocity profile determined at WUSU station also did not contain a midcrust low velocity zone. With these data we were unable to discern the eastern extent of this low velocity region; however, receiver functions sampling the crust to the south of DRAG are not indicative of a midcrust low velocity anomaly. The velocities and extent of the low velocity anomaly can be further constrained using the analysis scheme of Julia´ et al. (2000), which combines surface wave dispersion with the analysis of receiver functions. 3.3. Crustal Poisson’s ratio The v p/v s values obtained in the slant stacking analysis yield Poisson’s ratios between 0.20 and 0.30. While the errors in the Poisson’s ratios are of the order of this range, the values are fairly consistent for similar locations (Fig. 7). The results suggest that the bulk crustal Poisson’s ratio is slightly lower under the
173
Tianchi and Wangtian volcanoes than in regions immediately surrounding the volcanic craters, however, with these data, the trend is poorly constrained. The absolute values of v p/v s are not constrained by the NA direct search method. However in most cases, the change in v p/v s from the lower to upper crust is pronounced, consistent among multiple profiles determined at the same station, and had some spatial consistency: stations to the west of the volcanic region (JIYU, PSIA, JILA, LIAO) had a transition from a relatively low v p/v s in the upper crust to a higher v p/v s in the lower crust, while the opposite trend was resolved to the south and north of the volcanic region (DAND, DRAG, DHAA, FROG, JANA, THAA, and YABA; Fig. 14). At BACH and CANY stations, where we resolved the midcrust low velocity layer, the v p/v s profiles have a differing character (Fig. 14). The v p/v s ratio increases from upper to lower crust at CANY, while at BACH it decreases. With only two models, it is difficult to resolve this discrepancy, and is discussed further below.
4. Conclusions In this paper, we presented the crustal structure in the Changbaishan region based on modeling teleseismic receiver functions. The Moho depths we obtained in this study are consistent with those obtained by other researchers (Shin and Baag, 2000; Zhang et al., 2002; Song et al., in preparation). Researchers modeling wide-angle reflection/refraction data also detected a relatively thick crust (38–40 km) and a low velocity layer in the midcrust under the center of the Changbaishan volcanic region (Zhang et al., 2002; Song et al., in preparation). The low velocity anomaly extends laterally from about 30 km south of the Tianchi crater to about 50 km to the north, while the top of the anomaly is about 10 km below the surface and is about 10 km thick (Zhang et al., 2002; Song et al., in preparation). Furthermore, a magnetotelluric profile centered on the Tianchi crater and extending approximately north–south reveals a low resistivity anomaly extending from under the Tianchi to about 30 km to the north, at a depth of slightly over 10 km (Tang, 1998). Both of these results generally correspond to the low
174
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
velocity anomaly we resolved in this study. The fact that the region of low resistivity does not extend south of the Tianchi crater may provide a resolution between the contrasting v p/v s profiles determined for the crust south of BACH and of CANY, and may be related to an increased water content in the uppercrust between BACH station and the Tianchi compared to the crust under the Wangtian volcano. Indeed, most of the geothermal springs in this region are located to the north of the Tianchi (e.g., Shangguan and Sun, 1997). Based on these results, it is difficult to determine what the low velocity anomaly indicates. A likely explanation is that the decreased seismic velocities indicate elevated temperatures (e.g., Fountain and Christensen, 1989). The v p/v s results from BACH indicate that the low velocity anomaly may have an associated high Poisson’s ratio, which might indicate an increase in crack density and fluid content (e.g., O’Connell and Budiansky, 1974), possibly partial melt at depth (e.g., Mavko, 1980). Wilkens et al. (1984) and Brantley et al. (1990) showed that below 3.7–7 km depth the v p/v s changes little due to crack density, thus the increase in v p/v s at the depths we resolved may indicate partial melt. There is geothermal activity here, indicating elevated temperatures at depth, which could result in the presence of partial melt. Owens and Zandt (1997), using a receiver function analysis, detected Poisson’s ratios of 0.31– 0.33 in the lower crust of northern Tibet and argued for the presence of partial melt. We did not find Poisson’s ratios nearly as high in this study, although, we could not constrain the absolute values of Poisson’s ratio. Increasing mafic content tends to elevate v p/v s as well (e.g., Christensen and Fountain, 1975), and increased hydrothermal activity in the crust sampled by receiver functions recorded at BACH stations would also elevate v p/v s (O’Connell and Budiansky, 1974). With only these data we can not definitely resolve the presence of partial melt; however, we surmise that it is not present in large amounts. Finally, recent InSAR (interferometric synthetic aperture radar) observations in the Changbaishan region indicate a broad region of rapid subsidence (possibly up to 9 cm/year) essentially coincident with the lateral extent of the low velocity anomaly inferred in this study (Kim et al., 2001). This may imply that the hot region in the midcrust is
cooling rapidly due to circulation of meteoric water in a broad region of heavily cracked rock.
Acknowledgments This was a cooperative project between the Research Center of Geophysical Exploration, Chinese Seismological Bureau, and the State University of New York at Binghamton, covered under the Sino-US Seismological Research Protocol signed between USGS/NSF and CSB in 1980. We would like to thank IRIS/PASSCAL for providing the seismometers, as well as software, technical assistance and field support. Installation of the broadband network would not have been possible without the field assistance of Noel Barstow of PASSCAL and IRIS intern Margaret Boettcher. PASSCAL software assistance was provided by Sid Hellman and Paul Friberg. We thank Charles Ammon and Malcolm Sambridge for providing software used in this study. Most figures and maps were produced using GMT (Wessel and Smith, 1991). Careful reviews by F. Darbyshire, L. Pujades and editor K. Furlong are greatly appreciated, and EAH thanks R. van der Hilst and B. Hager for assistance and advice. The coordination and cooperation of the Changchun Office of the CSB and the staff of the RCGE, under the leadership of Dr. Zhang, was essential for the installation, operation and success of the temporary network. This project was funded under DOE grant DE-FC04-98A178918.
References Ammon, C.J., 1991. The isolation of receiver effects from teleseismic P waveforms. Bull. Seismol. Soc. Am. 81, 2504 – 2510. Ammon, C.J., Randall, G.E., Zandt, G., 1990. On the nonuniqueness of receiver function inversions. J. Geophys. Res. 95, 15303 – 15318. Brantley, S.L., Evans, B., Hickman, S.H., Crerar, D.A., 1990. Healing of microcracks in quartz: implications for fluid flow. Geology 18, 136 – 139. Chevrot, S., van der Hilst, R.D., 2000. The Poisson ratio of the Australian crust: geological and geophysical implications. Earth Planet. Sci. Lett. 183, 121 – 132. Christensen, N.I., 1996. Poisson’s ratio and crustal seismology. J. Geophys. Res. 101, 3139 – 3156. Christensen, N.I., Fountain, D.M., 1975. Constitution of the lower continental crust based on experimental studies of
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175 seismic velocities in granulite. Geol. Soc. Amer. Bull. 86, 227 – 236. Cui, Z., Wei, H., Liu, R., 1995. A textural research of record history about eruption activity of the Tianchi volcano. Volcanic Active and the Human Environment. Seismological Press, Beijing, pp. 36 – 39 (in Chinese). Darbyshire, F.A., Priestley, K.F., White, R.S., Stefa´nsson, R., Gudmundsson, G.B., Jakobsdo´ttir, S.S., 2000. Crustal structure of central and northern Iceland from analysis of teleseismic receiver functions. Geophys. J. Int. 143, 163 – 184. Du, Z.J., Foulger, G.R., 2001. Variations in crustal structure across central Iceland. Geophys. J. Int. 145, 246 – 264. Fan, Q., Liu, R., Li, D., Li, Q., 1999. Significance of K–Ar age of bimodal volcanic rocks at Wangtian volcano, Changbaishan Area. Chin. Sci. Bull. 44, 660 – 664 (in Chinese). Fountain, D.M., Christensen, N.I., 1989. Composition of the continental crust and upper mantle; a review. Geol. Soc. Amer. Bull. 172, 711 – 742. Guo, L.C., Ma, S.Z., Zhang, Y.S., 1996. The magma study of Changbai volcano by CT technique. CT Theory Appl. 5, 47 – 52 (in Chinese). Jin, B., Wei, X., 1994. Researching Volcanic Geology in Changbai Mountain. Northeast Korean Nationality Education Press, Yanbian (in Chinese). Julia`, J., Vila, J., Macia`, R., 1998. The receiver structure beneath the Ebro Basin, Iberian Peninsula. Bull. Seismol. Soc. Am. 88, 1538 – 1547. 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. Kennett, B.L.N., Engdahl, E.R., 1991. Travel times for global earthquake location and phase identification. Geophys. J. Int. 105, 429 – 465. Kim, S.W., Won, J.S., Kim, J.W., Moon, W.M., Min, K.D., 2001. Multi temporal JERS-1 SAR investigation of Mt. Baekdu stratovolcano using differential interferometry. Geosci. J. (Seoul) 5, 301 – 312. Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J. Geophys. Res. 84, 4749 – 4762. Liu, J., 1983. Study on Cenozoic Volcanic Activity in Changbaishan Area. Beijing Science and Technology Press, pp. 343 (in Chinese). Liu, R., Chiu, S.H., Cai, L.Z., Wei, H.Q., Yang, Q.F., Xi, Z.Q., Fu, G.C., Zhong, J., 1997. Chronological study and its implications of recent major eruptions of Tianchi volcano in Changbai Mountain. Sci. China 27, 437 – 441 (in Chinese). Liu, R.X., Wei, H.Q., Lee, J.T., 1998. The Recent Eruptions of the Tianchi Volcano in the Changbaishan. Science Publishers, Beijing, China. (in Chinese). Mavko, G.M., 1980. Velocity and attenuation in partially molten rocks. Geophys. J. Int. 85, 5173 – 5189. O’Connell, R.J., Budiansky, B., 1974. Seismic velocities in dry and saturated cracked solids. J. Geophys. Res. 79, 5412 – 5426.
175
Owens, T.J., Zandt, G., 1997. Implications of crustal property variations for models of Tibetan plateau evolution. Nature 387, 37 – 43. Phinney, R.A., 1964. Structure of the Earth’s crust from spectral behaviour of long-period body waves. J. Geophys. Res. 69, 2997 – 3017. Sambridge, M., 1999. Geophysical inversion with a neighborhood algorithm: i. Searching a parameter space. Geophys. J. Int. 138, 479 – 494. Shangguan, Z., Sun, M., 1997. Mantle-derived rare-gas releasing features at the Tianchi volcanic area, Changbaishan mountains. Chin. Sci. Bull. 42, 768 – 771. Shangguan, Z., Zheng, Y., Dong, J., 1997. Material sources of escaped gases from the Tianchi volcanic geothermal area, Changbai Mountains. Sci. China 40, 390 – 397. Sheriff, R.E., Geldart, L.P., 1982. Exploration Seismology: Volume 1. History, Theory and Data Acquisition. Cambridge University Press, Cambridge. Shin, J.S., Baag, C.E., 2000. Moho depths in the border region between northern Korea and northeastern China from waveform analysis of teleseismic pMP and pP phases. Geosci. J. (Seoul) 4, 313 – 320. Song, J.L., Wu, F.T., Hetland, E.A., Zhang, X.K., 2003. Wide angle refraction/reflection P wave velocity structure of Changbai Volcano, China. In preparation. Tang, J., 1998. The electronic structure of Tianchi volcano area by magnetotelluric. In: Liu, R.X., Wei, H.Q., Lee, J.T. (Eds.), The Recent Eruptions of the Tianchi Volcano in the Changbaishan. Science Publishers, Beijing, China, pp. 108 – 123 (in Chinese). van der Hilst, R.D., Kennett, B.L.N., Shibutani, T., 1998. Upper mantle structure beneath Australia from portable array deployments. In: Braun, J. (Ed.), Structure and Evolution of the Australian Continent, volume 26 of Geodynamics. American Geophysical Union, Washington, DC, pp. 39 – 57. Wessel, P., Smith, W.H.F., 1991. Free software helps map and display data. EOS Transactions, American Geophysical Union 72, 441, 445–446. Wilkens, R., Simmons, G., Caruso, L., 1984. The ratio of vp/vs as a discriminant of composition for siliceous limestones. Geophysics 49, 1850 – 1860. Zhang, J., Langston, C.A., 1995. Dipping structure under Doubes, Belgium, determined by receiver function modeling. Bull. Seismol. Soc. Am. 85, 254 – 268. Zhang, C.K., Zhang, X.K., Zhao, J.R., Liu, B.F., Zhang, J.S., Yang, Z.X., Hai, Y., Sun, G.W., 2002. Crust–mantle structure of the Changbaishan Tianchi volcanic region and its vicinity: an exploratory study and inferences. Chin. J. Geophys. 45, 862 – 871. Zhu, L., Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver functions. J. Geophys. Res. 105, 2969 – 2980.
Crustal structure in the Changbaishan volcanic area, China, determined by modeling receiver functions E.A. Hetlanda,*, F.T. Wub, J.L Songb,1 a
Department of Earth, Atmospheric and Planetary Sciences, Massachusetts Institute of Technology, 77 Massachusetts Avenue, Cambridge, MA 02139-4307, USA b Department of Geological Sciences and Environmental Studies, State University of New York at Binghamton, Binghamton, NY 13902, USA Received 21 July 2003; accepted 2 June 2004 Available online 22 July 2004
Abstract During 1998–1999, we installed a temporary broadband seismic network in the Changbaishan volcanic region, NE China. We estimated crustal structure using teleseismic seismograms collected at the network. We detected a near surface region of strong anisotropy directly under the main volcanic edifice of the volcanic area. We modeled 109 receiver functions from 19 broadband stations using three techniques. First we used a bslant-stackingQ method to model the principal crustal P reverberation phases to estimate crustal thickness and the average crustal P to S speed ratio (v p/v s), assuming an average P-wave velocity in the crust. We then estimated crustal S-wave velocity (v s) and v p/v s profiles by modeling stacked receiver functions using a direct search. Finally, we inverted several receiver functions recorded at stations closest to the main volcanic edifice using least squares to estimate v s velocity profiles, assuming a v p/v s value. The results from the three estimation techniques were consistent, and generally we found that the receiver functions constrained estimates of changes in wave speeds better than absolute values. We resolved that the crust is 30–39 km thick under the volcanic region and 28–32 km thick away from the volcanic region, with a midcrust velocity transition at about 10–15 km depth. We estimated that the average crust P-wave velocity is about 6.0–6.2 km/s surrounding the main volcanic region, while it is slightly lower in the vicinity of the main volcanic edifice. The estimates of v p/v s were more ambiguous, but we inferred that the bulk crustal Poisson’s ratio (which is related to v p/v s) ranges between 0.20 and 0.30, with a suggestion that the Poisson’s ratio is lower under the central volcanic region compared to the surrounding areas. We resolved low S-wave velocities (down to about 3 km/s) in the middle crust in the region of the main volcanic edifice. The low velocity anomaly extends from about 5–10 to 15–25 km below the surface, probably indicating a region of elevated temperatures. We were unable to determine if partial melt is present with the data we considered in this paper. D 2004 Elsevier B.V. All rights reserved. Keywords: Changbai Mountains; PASSCAL; Crustal structure; Receiver functions; Tianchi volcano
* Corresponding author. E-mail address: [email protected] (E.A. Hetland). 1 Now at US Geological Survey, 384 Woods Hole Road, Woods Hole, MA 02543 USA. 0040-1951/$ - see front matter D 2004 Elsevier B.V. All rights reserved. doi:10.1016/j.tecto.2004.06.001
158
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
1. Introduction The Changbaishan volcanic region is located in Jilin Province in northeastern China, along the border between China and North Korea (Fig. 1). There are two main volcanic edifices in the Changbaishan region, the Wangtian and the Tianchi, which are of basaltic and andesitic to rhyolitic composition underlain by a large Cenozoic basaltic shield (Liu, 1983; Liu et al., 1998). The Wangtian was last active 2.12 Ma and has approximately 400 km2 of exposed volcanic rock associated with it (Fan et al., 1999). The Tianchi volcano became active about 1.66–1.48 Ma years ago (Jin and Wei, 1994) and is still active, with the most recent activity in 1903 (Cui et al., 1995). The Tianchi eruption of 1000AD was one of the largest recorded in China in the past 2000 years (Liu et al., 1997). Until recently, little geophysical research has been conducted in the Changbaishan region. Regional Swave tomography by Guo et al., (1996) indicated a low velocity anomaly at 25–75 km depth below the Tianchi volcano (also see Liu et al., 1998). Additionally, a magnetotelluric study done by Tang (1998) revealed a region of low resistivity below 10–15 km depth extending from the Tianchi volcanic crater to about 30 km to the north. The region of low resistivity might indicate a region of partial melt or a region of high fluid content, but was interpreted to be most probably a magma chamber by Tang (1998). The Moho depth was determined to be about 30–35 km by modeling surface and Moho reflected phases from a regional earthquake (Shin and Baag, 2000). Based on geochemical analysis, Shangguan and Sun (1997) have determined that the gases released from local hot springs are from the top of a residual, mantle-derived magma chamber in the crust, located possibly as shallow as 5 km below the Tianchi crater. The water from hot springs in the Changbaishan region are about 80.75 8C and have been determined to be meteoric water circulating through regions of hot rock (Shangguan et al., 1997). In the summer of 1998, the State University of New York at Binghamton, in collaboration with the Research Center of Geophysical Exploration (RCGE) of the Chinese Seismological Bureau, installed 19 portable, broadband seismic stations in NE China, in
the region of the Changbaishan volcanic area (Table 1 and Fig. 1). The stations were borrowed from IRIS PASSCAL and each station comprised a Guralp 3T sensor and a six-channel Ref Tek recorder. We established stations at 18 sites in late June 1998. Eight sites were occupied until late September 1998, while 10 sites were occupied until March 1999. The one remaining station was installed in September 1998 and removed in March 1999. In August 1998, the RCGE conducted an active source, wide-angle reflection/refraction experiment, and analysis of profiles near the Tianchi volcano indicated a low Pwave velocity anomaly under the Tianchi volcano at about 10–20 km depth, and resolved fairly uniform Moho depths of about 30–35 km away from the volcanic center, thickening to about 40 km under the Tianchi edifice (Zhang et al., 2002; Song et al., in preparation). In this study, we present a determination of crustal structure using teleseismic receiver functions throughout the Changbaishan region. The receiver function method has been applied to many tectonic problems (e.g., Owens and Zandt, 1997; Julia` et al., 1998; van der Hilst et al., 1998; Chevrot and van der Hilst, 2000; Julia´ et al., 2000; Zhu and Kanamori, 2000) and volcanic regions (e.g., Langston, 1979; Darbyshire et al., 2000; Du and Foulger, 2001).
2. Receiver functions 2.1. Source equalization The first 20 or so seconds following direct P in a teleseismic seismogram principally contains rupture effects, source-side structure, and receiver-side structure. The source-side signal results from near source reflected P waves and S to P conversions. Receiverside signal consists of P to S conversions and P and S reflections near the receiver, principally in the crust (Fig. 2). Direct P is incident on the Moho nearly vertically for a teleseismic distance source. Therefore, the amplitude of S conversions are larger on the horizontal components, while the P signals are larger on the vertical component. Since the horizontal components of the seismogram principally contain receiver-side effects, it is possible to equalize the source effects and retain the receiver-side effects by
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
159
Fig. 1. Map of the temporary broadband station locations, circle type indicates the approximate deployment of each station. The Tianchi and Wangtian volcanic craters are indicated by white triangles, rivers, lakes and ocean are black, and the border between China and North Korea is a thick dashed line.
deconvolving the vertical component from the radial component leaving a radial receiver function (e.g., Phinney, 1964; Langston, 1979; Ammon, 1991). In
order to avoid division by small numbers, the equalization procedure of Langston (1979), modified by Ammon (1991), uses a water-level deconvolution
160
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
Table 1 Locations of the temporary broadband stations, the geographic region of each station relative to the Tianchi volcano (bTianchiQ signifies the station is near the main volcanic edifice; see Fig. 13), and the number of receiver functions (No.) determined at each station used in this study Station
Latitude (8N)
Longitude (8E)
Region
No.
Station
Latitude (8N)
Longitude (8N)
Region
No.
BACH CBAI DAND DRAG JANA JIYU LING PSIA TILL YABA
42.514 42.058 40.146 42.504 41.362 42.379 41.799 42.950 42.267 43.003
128.158 128.063 124.391 129.031 126.512 126.795 126.957 126.078 123.914 129.489
Tianchi Tianchi west north south south south west west north
5 1 8 5 19 5 3 8 3 7
CANY CH2A DHAA FROG JILA LIAO MANG THAA WUSU
41.951 43.801 43.344 41.420 43.728 41.272 41.948 41.697 42.326
128.005 125.450 128.198 128.174 126.664 123.284 127.594 125.998 127.277
Tianchi west north south west west south south south
3 3 16 7 3 5 3 3 2
where signal with spectral power below some threshold, called the water level, is discarded. This equalization method includes a Gaussian filter in the deconvolution, which acts as a low-pass filter removing noise in the seismograms. Because the receiver function method uses crustal converted and reflected phases to determine crustal structure, the region of crust sampled depends on the offset of these crustal reverberations (Fig. 2). The offsets of the reverberations depend on the ray parameter ( p) and the crustal thickness (or more generally the depth to the deepest interface). We estimated p from the event depth and range using the IASP91 model of P global seismic travel times (Kennett and Engdahl, 1991). The lateral extent from the station that the crust is sampled by a receiver function is at least the distance to the penetration point of the first bounce P arrival (Fig. 2). In this region where the crust is expected to be about 30–35 km thick, with an average P wave velocity of 6.0–6.5 km/s, and for the data set we used, the receiver
functions sample the Moho about 40–50 km from the stations. Furthermore, Sheriff and Geldart (1982) have demonstrated that the lateral resolution of seismically imaged interfaces is approximated by the radius of the Fresnel zone, and as the radius of the Fresnel zone increases, the P wave samples a larger area of an interface. Since the radius of the Fresnel zone increases with depth, the lateral resolution over an interface decreases with depth. In this area, the lateral resolution of P on the Moho is about 10–15 km. Energy propagating off the back azimuth will decrease the amplitude of the SV waves on the radial component and will not be resolved on the equalized radial receiver function (e.g., Ammon, 1991; Zhang and Langston, 1995). For this reason it is crucial to rotate the seismogram accurately, and strong lateral heterogeneity and receiver-side anisotropy will complicate this. In practice, equalization is not perfect, thus we considered all the equalized receiver functions to be possibly corrupted with noise, and we only
Fig. 2. Cartoon showing direct P (labeled Pp) and the dominant crustal reverberations, where, except for mantle P, lowercase and uppercase indicate upgoing and downgoing rays, respectively (a), and the associate receiver function (b); after Ammon et al. (1990).
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
tried to model dominant arrivals and arrivals which were consistent for all traces.
161
band seismometers were located south of the network, since the highest concentration of seismicity meeting the above criterion was in the seismic belt from Indonesia through Papua New Guinea to the Fiji islands (Fig. 3). We used one event from the Mariana Trench. We found only three events from the north (the Aleutian Trench region) and one from the west (located in Iran). Due to equipment failures, power outages, periods of unusually high anthropogenic noise levels, and differing durations of deployment not all of the events were recorded at all of the stations. In this study, we used 109 seismograms recorded at the 19 stations (Table 1). We verified the orientations of all stations, as well as the back azimuths of all events, by rotating the unequalized horizontal seismograms and visually inspecting the P energy on the radial and tangential
2.2. Data In this study, out of 41 earthquakes considered, we used 29 events below 48 km depth (one event was shallow, with a depth listed as 33 km) and at a range of 30–908 (Table 2 and Fig. 3). The first criterion was used to avoid complex near-source structure effects, while the second was used to ensure that the direct P arrived well before secondary P phases. We also used only large events to ensure that the signal was strong. The minimum magnitude event we found to equalize well was 5.3 Mb, although this depended on the signal-to-noise ratio of a given seismogram. Most of the accepted events detected by the temporary broad-
Table 2 Earthquakes used in this study Year
Day
Time
Latitude (8N)
Longitude (8E)
Depth (km)
Magnitude
Location
1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1998 1999 1999 1999 1999 1999 1999 1999 1999 1999 1999
190 190 197 206 231 240 245 245 248 251 258 263 264 271 284 300 319 329 361 019 028 028 035 036 037 039 061 063 075
14:45:39.9 19:39:43.9 11:56:36.4 02:39:23.3 17:41:43.6 12:40:58.7 08:37:29.9 18:52:42.2 17:09:31.4 09:10:03.0 08:35:51.5 21:21:55.1 06:52:41.1 13.34:30.4 12:04:54.7 21:16:21.0 02:44:12.3 18:05:25.7 00:38:26.7 03:35:33.8 08:10:05.4 18:24:25.2 19:28:00.8 11:39:45.1 21:47:59.4 16:22:54.3 21:15:30.0 05:38:26.5 03:32:16.2
30.487 60.530 11.040 13.608 6.179 0.154 5.410 29.695 4.051 13.257 5.624 7.770 0.262 8.194 21.040 2.921 21.589 7.859 21.632 4.596 52.886 4.582 4.033 12.616 12.853 7.908 51.870 28.343 2.683
178.994 153.219 166.160 166.867 147.745 125.018 126.764 178.790 152.041 144.007 151.637 106.951 122.467 112.413 179.110 128.623 176.504 158.622 176.376 153.235 169.123 153.658 95.283 166.966 166.697 147.734 179.430 57.193 125.674
130.0 145.0 110.0 44.0 55.0 66.0 50.0 230.0 192.0 141.0 83.0 67.0 147.0 152.0 624.0 61.0 149.0 48.0 144.0 114.0 67.0 101.0 56.0 213.0 90.0 70.0 72.2 33.0 114.0
6.4ML 6.2ML 6.9ML 6.2ML 5.3MB 6.2MB 6.8ML 5.4MB 5.5MB 5.8MB 6.2ML 5.5MB 6.1MB 6.5ML 5.9ML 5.9MB 6.2ML 6.2ML 6.6ML 6.5ML 6.6ML 6.3ML 5.9ML 5.9ML 7.3MS 5.9ML 5.5MB 6.5MS 5.6MB
KERMADEC ISLANDS, NEW ZEALAND SOUTHERN ALASKA SANTA CRUZ ISLANDS VAUATU ISLANDS EASTERN NEW GUINEA REG., P.N.G. SOUTHERN MOLUCCA SEA MINDANAO, PHILIPPINE ISLANDS KERMADEC ISLANDS, NEW ZEALAND NEW BRITAIN REGION, P.N.G. MARIANA ISLANDS NEW BRITAIN REGION, P.N.G. JAWA, INDONESIA MINAHASA PENINSULA, SULAWESI JAWA, INDONESIA FIJI ISLANDS REGION HALMAHERA, INDONESIA FIJI ISLANDS REGION SOLOMON ISLANDS FIJI ISLANDS REGION NEW IRELAND REGION, P.N.G. FOX ISLANDS, ALEUTIAN ISLANDS NRE IRELAND REGION, P.N.G. NORTHERN SUMATERIA, INDONESIA SANTA CRUZ ISLANDS SANTA CRUZ ISLANDS EASTERN NEW GUINEA REG., P.N.G. A RAT ISLANDS, ALEUTIAN ISLANDS SOUTHERN IRAN TALAUD ISLANDS, INDONESIA
Note that the Celebes Sea region includes the Philippine Islands and Molucca Sea events.
162
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
Fig. 3. Epicentral locations (circles) of the earthquakes used in this study and great circle path (dashed lines) from the epicenter to the seismic array (approximated as the location of the Tianchi volcano). The square indicates the area of Fig. 1, and the black lines are major plate boundaries.
components. For a typical seismogram, when rotated correctly, all of the energy of direct P is on the radial component until about 5 s later when the Ps phase (S conversion at the Moho) is detected (Fig. 4a). We were unable to rotate all of the energy off the tangential component of seismograms recorded at CBAI, where there is equal amplitude energy on both the tangential and radial components in the first several seconds (Fig. 4b). Moreover, this energy appears to be out of phase on the two components, which is typical of anisotropy observations. Only seismograms of events to the south of CBAI were affected, and this led us to conclude that there is a region of possible anisotropy in the crust to the south of CBAI. We observed only a small hint of this phenomenon on seismograms from CANY, located about 10 km to the south of CBAI, and since the seismic waves to both of these stations sampled
similar lower crust but more distinct upper crust, we surmised that the region of anisotropy is in the shallow crust under the Tianchi crater. The anisotropy may be the result of cracked rock corresponding to the geothermal activity here, and it may be coincident with the source of gases detected in hot springs in this area. Due to this complexity in the seismograms, we only used one seismogram recorded at CBAI, which was the record of an event, located to the north, that was not affected by the anisotropy. We do not present an analysis of the region of potential anisotropy in this paper. 2.3. Analysis After identifying all usable seismograms, we determined the receiver functions using a constant Gaussian width parameter of 2.4 Hz and varying
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
163
Fig. 4. Radial and tangential components of seismograms of two events located in the Santa Cruz islands recorded at two stations: (a) FROG, located about 70 km to the south of the Tianchi crater and (b) CBAI, located adjacent and slightly to the north of the Tianchi. Time is relative to the approximate arrival time of the initial P wave, and events are to the south of the stations.
water levels (from 510-5 to 10-2). We determined the water level largely by the quality of the original seismogram and stability of the equalization. We interpreted the radial receiver functions using three analysis techniques which estimate differing sets of model parameters: a bslant-stackingQ procedure, a direct search forward modeling scheme, and a leastsquares inversion (LS). In all of these methods, we assumed that the crust is locally homogeneous, isotropic, and that the Moho is not dipping. Since the receiver functions sample relatively localized regions of the crust, we only assumed the crust to be homogeneous over these scales. We could not account for anisotropy in these analyses, but we did not use any seismograms that showed signs of
anisotropy. We used the results of the three methods to cross-check each other. While using the above three techniques does not provide independent measurements, due to the differing methods and model parameters estimated, the three methods indicate robustness of the final results. 2.3.1. Slant-Stacking The slant-stacking method has the ability to resolve simultaneously crustal thickness (H) and the P to S wave speed ratio (v p/v s), thereby partially resolving the thickness/velocity tradeoff (Zhu and Kanamori, 2000; Chevrot and van der Hilst, 2000). In this calculation, only three travel times, relative to direct P, are needed: Ps (P–S conversion at the Moho, where,
164
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
except for mantle P, lowercase and uppercase indicate upgoing and downgoing rays, respectively, in the crust; Fig. 2), PpPs, and PpSs (or PsPs). Following the
arguments of Chevrot and van der Hilst (2000), we weighted each of these arrivals the same. Chevrot and van der Hilst (2000) did not use the PpSs phase, but
Fig. 5. (a) Receiver functions determined at CANY station, with the predicted travel time curves estimated from the optimum solution of H and v p/v s determined by slant-stack analysis (v¯ p=6.0 km/s). (b) Solution surface constructed by a direct search over H and v p/v s with the optimum solution (star) and the 0.9 contour line. Scale is normalized and all solutions with negative sum are white.
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
we found including it better constrained the solution with these data. The differential travel-time equations depend on H, v p/v s, the average crustal P-wave
165
velocity (v¯ p) and the ray parameter ( p). The ray parameter is determined from the event depth and range, and we assumed values of v¯ p of 6.0 and 6.5 km/s.
Fig. 6. (a) Receiver functions determined at JANA station, with the predicted travel time curves estimated from the optimum solution of H and v p/v s determined by slant-stack analysis (v¯ p=6.5 km/s). (b) Solution surface constructed by a direct search over H and v p/v s with the optimum solution (star) and the 0.9 contour line. Scale is normalized and all solutions with negative sum are white.
166
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
Poisson’s ratio (m) can be determined from v p/v s through the relation ! 1 1 1 m¼ 2 2 vp =vs 1 and is a useful quantity for crustal studies since m is sensitive to mineral composition (e.g., Christensen, 1996), the presence of cracks and fluid content (e.g., O’Connell and Budiansky, 1974), and partial melt (e.g., Mavko, 1980). At each station, for a given choice of H and v p/v s, we summed the amplitudes of the receiver functions, at the predicted differential travel times of the Ps, PpPs, and PpSs phases: X s H; vp =vs ¼ aðtPPs Þ þ a tPPpPs a tPPpSs R
where R is a set of receiver functions and a(t) is the amplitude of a receiver function at time t. We grid searched over (H,v p/v s) parameter space, and we took the optimum solution to be at the maximum of the solution surface (Figs. 5 and 6). In the case of a
sufficient number of receiver functions and clear arrivals, the model space will show one global maximum (Fig. 6). In the neighborhood of the maximum, the parameter surface (s) has a trend which is related to the v p/v s–H tradeoff. Previously Zhu and Kanamori (2000) used the curvature of s in the two parameter directions to infer errors on the estimated parameters. Noting that this method does not account for the anticorrelation between v p/v s and H, Chevrot and van der Hilst (2000) used a boot-strapping method, where they systematically removed several traces and determined the solution for each realization, then taking the standard deviation of these solutions to estimate the errors. While this method of error estimation accounts for the obvious anticorrelation, we found that in this study it underestimated the errors, possibly due to the relatively small number of traces we had for any given station, the data quality, or the equalization method we used. We determined that the boot-strapping method underestimated the errors through forward trials, where we perturbed the optimum solution until the change in travel-time curves was beyond the resolution of the data. We
Fig. 7. Estimated crustal thickness (H, top) and Poisson’s ratio (v, bottom) from the slant-stacking analysis, using v¯ p=6.0 and 6.5 km/s, as a function of distance from the Tianchi crater. Error bars are the minimum and maximum values of H or v of the error ellipses (e.g., Figs. 5 and 6).
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
found that the 90% contour line of s was an effective proxy for the errors (Figs. 5 and 6). Since information on v¯ p is scant in this area, we performed slant stacking using both 6.0 and 6.5 km/s. While the values of H and m changed according to the value of v¯ p we used, the trends of changes in H and m across the broadband network remained consistent and the values were consistent with neighboring stations (Fig. 7). Due to the small number of receiver functions determined at CBAI and CH2A, we did not include these stations in the slant-stacking analysis. 2.3.2. Direct search We also modeled the observed receiver functions with a direct search algorithm, which uses a nearest-
167
neighbor search method (Sambridge, 1999). The nearest-neighbor algorithm directs a Monte Carlo search to a global solution, and in that it is similar to a genetic algorithm or a simulated annealing algorithm. The nearest-neighbor algorithm (NA) has been shown to be effective for complex problems and has been applied successfully to receiver function modeling, solving for v p/v s, layer thicknesses and v s at the top and bottom of each model layer (Sambridge, 1999). As opposed to other direct search algorithms, the NA is only dependent on one tuning parameter. Due to complexity of search space and number of parameters per layer, in practice we were unable to get good results using more than four or at most five model layers. Thus, the NA tended to pick up the
Fig. 8. Stack of all of the receiver functions determined at CANY station from events in the Celebes Sea region (thin line, top panel) with the best fit model (dashed line). The best fit model (solid line, bottom panels) was found using the NA direct search. Fit is normalized such that the best fitting model has misfit of 1.0, and the all of the models with at least twice the misfit of the best model are shown (white to gray lines) along with the average (dashed line) of the shown models.
168
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
major velocity discontinuities, while resolving velocity gradients between the discontinuities (Fig. 8). Solving for v p/v s is a large asset to this method. We found the NA to be more sensitive to changes in v p/ v s, which indicate horizons of changing properties, and less sensitive to the actual value of v p/v s (e.g., Fig. 9). We inverted stacked receiver functions using the NA. By stacking several receiver functions with similar back azimuths and ray parameters, it was possible to increase the signal-to-noise ratio of receiver functions. Due to the small number of receiver functions determined at most stations, we could not get a stacked receiver function for LING, MANG, TILL, or WUSU. At TILL station we modeled a single receiver function which had a particularly high signal-to-noise ratio. At most of the other stations there were enough receiver func-
tions to get at least one stacked function to use in this analysis, which we modeled to determine velocity profiles (Fig. 8). We determined three velocity profiles from three stations on the outskirts of the volcanic region (BRIT, DHAA, and JANA; Fig. 9), and two velocity profiles for stations DAND, JIYU, and YABA. Except for the event used at TILL, all events we used in this phase of the study were from south to southeast directions. So that resulting velocity profiles did not contain a discontinuity when the data did not require one, we searched over a large range of layer thicknesses and specified the same range of velocities at the top and bottom of adjacent layers. While we were unable to model the entire receiver function using a model with 4–5 layers, we were able to model the dominant arrivals (Figs. 8 and 9). In all cases, the Moho depths determined with the NA direct search were consistent
Fig. 9. Stacked receiver functions modeled using a NA direct search. Receiver functions were determined at JANA station and are from the New Britain Islands, Celebes Sea and Fiji Islands regions. Caption is as in Fig. 8.
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
169
Fig. 10. Comparison of the crustal thickness (H) and average v p (v¯ p) estimated at CANY and JANA stations by the slant-stack and NA direct search methods. Diamond and black crosses are the H and v¯ p determined by the NA method (diamond represents the one estimate at CANY, and size of cross indicates standard deviation of the three estimates at JANA). Gray circles are H estimated in the slant-stack analysis for the indicated v¯ p, error bars are shown as in Fig. 7.
with those obtained in slant-stacking, if the differences in average crustal v p were accounted for (Fig. 10). Additionally, at the stations where we were able to determine more than one velocity model, the results were consistent (Fig. 9).
2.3.3. Least squares The last technique we used to model the receiver functions was the least-squares inversion of Ammon et al. (1990). The least-squares inversion (LS) method requires an initial velocity model composed of an
Fig. 11. Two receiver functions (thin lines, left panels) determined at CANY station from events in the Celebes Sea region on days 240 (cele240) and 245 (cele245) of 1998, and the best fit modeled receiver functions determined through LS inversion (dashed thick lines, right panels) using a starting model composed of 40 one kilometer layers (initial velocity profile is thin gray line). The velocity models determined from modeling the stack of cele240 and cele245 with the NA direct search (Fig. 8) are shown for reference.
170
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
arbitrary number of layers. The v s in each layer is then estimated, holding the layer thicknesses constant and assuming a constant Poisson’s ratio of 0.25 (Ammon et al., 1990). We needed to apply smoothing constraints in the inversion in order to avoid checker-boarding in
final solution. The LS technique is more sensitive to velocity discontinuities than to the absolute magnitude of the velocities (Ammon et al., 1990). Our motivation for using this inverse technique was to determine higher resolution crustal velocity
Fig. 12. Velocity models determined by LS inversion of the receiver functions at stations closest to the Tianchi volcanic crater. Individual velocity models are shown as thin gray lines, and the average of all velocity profiles for each station is shown with thick lines. Gray lines projecting from the stations show the approximate extent from the stations that each receiver function sampled the crust along the great circle path. Cross-hatched regions show the approximate regions of the crust sampled by the receiver functions which detected a midcrust low-velocity anomaly. The locations of the Tianchi and Wangtian volcanic craters (triangles) and the Chinese national border (dashed line) are shown for reference.
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
profiles than the NA direct search allowed. We used a starting model with 40–45 1-km-thick layers. We applied varying levels of smoothing until we were able to model the dominant phases while avoiding
171
overcomplexity and checker-boarding in the resultant model. Fig. 11 shows the results of inverting two receiver functions determined from events in the Celebes Sea region recorded at CANY station, where
Fig. 13. Velocity models determined from modeling receiver functions recorded by the Changbaishan broadband network, grouped according to region. Average velocity models, across the indicated stations, determined using the NA direct search (smooth bold lines). The average velocity models from each NA direct search (solid, thin gray lines) are shown to indicate the spread. The average velocity models from the 40–45 layer LS inversion (stepped lines) are shown for reference. NA velocity models from BACH and CANY are the results of single waveform analysis. Gray lines projecting from the stations show the approximate extent from the stations that each receiver function sampled the crust along the great circle path, and right and left leaning striped backgrounds indicate v p/v s as in Fig. 14a and b, respectively.
172
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
3. Discussion
volcanic region, and 30–39 km in the volcanic region, with the thickest crust closest to the Tianchi crater (Figs. 7 and 13). All of the events recorded at CANY station are from the south to southeast, so the region of Moho sampled is to the south and southeast. Thus, receiver functions recorded at CANY station constrain crustal thickness south of the Tianchi crater, under the Wangtian volcanic edifice. Using direct search and least squares inversion results, we resolved the thickest crust in the area extending from FROG north to BACH, including possibly MANG (Figs. 12 and 13). We determined the v¯ p in the regions surrounding the main volcanic area (the Tianchi and Wangtian craters) to be about 6.0–6.2 km/s. The lowest v¯ p we found using the stacked waveforms were for BACH, CANY and YABA stations (all below 6.0 km/s with YABA being the lowest at v¯ pc5.7 km/s). Furthermore, all of the velocity profiles determined for stacked receiver functions contained a midcrust boundary at about 10–15 km depth, with the exception of YABA station, whose velocity profiles indicated a midcrust transition at 20–25 km depth (Fig. 13).
3.1. Crustal thickness and bulk velocity structure
3.2. Midcrust low velocity anomaly
Using the slant-stacking algorithm, we found the Moho to be at about 28–32 km depth away from the
The velocity profiles we obtained from modeling stacked receiver functions from BACH and CANY
the receiver function from the event on day 245 has a higher signal-to-noise ratio than the event on day 240. We only performed these higher resolution inversions using traces recorded at stations closest to the center of the volcanic region (BACH, CANY, DRAG, FROG, MANG, WUSU and one trace recorded at CBAI from the north, which was not as affected by anisotropy as the records from the south were). The crustal velocity structure of this region was of most interest and we wished to constrain the results of the NA direct search. We inverted a total of 26 traces recorded at these seven stations. In these inversions, we needed to apply relatively large smoothing (20– 40%) because of the overall data quality and noise level. There were five traces for which we could not obtain stable or consistent results and we discarded them. Furthermore, we only obtained one velocity profile at CBAI and WUSU stations. The resulting velocity profiles were consistent for similar areas of crust and are consistent with those determined above using the NA method (Fig. 12).
Fig. 14. Profiles of v p/v s determined through modeling stacked receiver functions through NA direct search. (a) Average v p/v s profiles (gray lines) determined at DAND, DRAG, DHAA, FROG, JANA, THAA and YABA stations (corresponding to the right leaning striped stations in Fig. 13) and the average of the set (thick line). (b) Average v p/v s profiles determined at LIAO, JILA, JIYU and PSIA stations (left leaning striped) and the average of the set (thick line). (c) Average v p/v s profiles determined from modeling stacked receiver functions at BACH and CANY.
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
both indicate a pronounced high velocity layer in the upper crust, with a low velocity anomaly in the middle crust (Fig. 12). The high velocity layer starts at about 5–10 km depth, and is about 10 km thick at BACH and 10–15 km thick at CANY. The S-wave velocities in the upper layer are over 4 km/s, while the velocities decrease to about 3 km/s or slightly lower in the middle crust. The high velocity layer is more pronounced in the velocity profile at CANY than at BACH, and the stacked receiver functions from both stations sampled crust to the south of the stations. We further investigated the crust in the volcanic region by modeling single receiver functions with least-squares inversions (Ammon et al., 1990). The velocity profiles indicated an ubiquitous low velocity zone in the midcrust which appears to be localized in the region extending from FROG to BACH station, possibly further to the north (Fig. 12). The low velocity anomaly is thicker and deeper south of CANY than under the north slope of the Tianchi volcano. Receiver functions that sampled crust near FROG station resolved a slight low velocity region at greater depth in the crust (Fig. 12). Receiver functions recorded at MANG do not suggest a low velocity anomaly and we concluded that the western extent of the low velocity anomaly is located under the western flank of the Changbaishan (Fig. 12). One velocity profile determined at WUSU station also did not contain a midcrust low velocity zone. With these data we were unable to discern the eastern extent of this low velocity region; however, receiver functions sampling the crust to the south of DRAG are not indicative of a midcrust low velocity anomaly. The velocities and extent of the low velocity anomaly can be further constrained using the analysis scheme of Julia´ et al. (2000), which combines surface wave dispersion with the analysis of receiver functions. 3.3. Crustal Poisson’s ratio The v p/v s values obtained in the slant stacking analysis yield Poisson’s ratios between 0.20 and 0.30. While the errors in the Poisson’s ratios are of the order of this range, the values are fairly consistent for similar locations (Fig. 7). The results suggest that the bulk crustal Poisson’s ratio is slightly lower under the
173
Tianchi and Wangtian volcanoes than in regions immediately surrounding the volcanic craters, however, with these data, the trend is poorly constrained. The absolute values of v p/v s are not constrained by the NA direct search method. However in most cases, the change in v p/v s from the lower to upper crust is pronounced, consistent among multiple profiles determined at the same station, and had some spatial consistency: stations to the west of the volcanic region (JIYU, PSIA, JILA, LIAO) had a transition from a relatively low v p/v s in the upper crust to a higher v p/v s in the lower crust, while the opposite trend was resolved to the south and north of the volcanic region (DAND, DRAG, DHAA, FROG, JANA, THAA, and YABA; Fig. 14). At BACH and CANY stations, where we resolved the midcrust low velocity layer, the v p/v s profiles have a differing character (Fig. 14). The v p/v s ratio increases from upper to lower crust at CANY, while at BACH it decreases. With only two models, it is difficult to resolve this discrepancy, and is discussed further below.
4. Conclusions In this paper, we presented the crustal structure in the Changbaishan region based on modeling teleseismic receiver functions. The Moho depths we obtained in this study are consistent with those obtained by other researchers (Shin and Baag, 2000; Zhang et al., 2002; Song et al., in preparation). Researchers modeling wide-angle reflection/refraction data also detected a relatively thick crust (38–40 km) and a low velocity layer in the midcrust under the center of the Changbaishan volcanic region (Zhang et al., 2002; Song et al., in preparation). The low velocity anomaly extends laterally from about 30 km south of the Tianchi crater to about 50 km to the north, while the top of the anomaly is about 10 km below the surface and is about 10 km thick (Zhang et al., 2002; Song et al., in preparation). Furthermore, a magnetotelluric profile centered on the Tianchi crater and extending approximately north–south reveals a low resistivity anomaly extending from under the Tianchi to about 30 km to the north, at a depth of slightly over 10 km (Tang, 1998). Both of these results generally correspond to the low
174
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175
velocity anomaly we resolved in this study. The fact that the region of low resistivity does not extend south of the Tianchi crater may provide a resolution between the contrasting v p/v s profiles determined for the crust south of BACH and of CANY, and may be related to an increased water content in the uppercrust between BACH station and the Tianchi compared to the crust under the Wangtian volcano. Indeed, most of the geothermal springs in this region are located to the north of the Tianchi (e.g., Shangguan and Sun, 1997). Based on these results, it is difficult to determine what the low velocity anomaly indicates. A likely explanation is that the decreased seismic velocities indicate elevated temperatures (e.g., Fountain and Christensen, 1989). The v p/v s results from BACH indicate that the low velocity anomaly may have an associated high Poisson’s ratio, which might indicate an increase in crack density and fluid content (e.g., O’Connell and Budiansky, 1974), possibly partial melt at depth (e.g., Mavko, 1980). Wilkens et al. (1984) and Brantley et al. (1990) showed that below 3.7–7 km depth the v p/v s changes little due to crack density, thus the increase in v p/v s at the depths we resolved may indicate partial melt. There is geothermal activity here, indicating elevated temperatures at depth, which could result in the presence of partial melt. Owens and Zandt (1997), using a receiver function analysis, detected Poisson’s ratios of 0.31– 0.33 in the lower crust of northern Tibet and argued for the presence of partial melt. We did not find Poisson’s ratios nearly as high in this study, although, we could not constrain the absolute values of Poisson’s ratio. Increasing mafic content tends to elevate v p/v s as well (e.g., Christensen and Fountain, 1975), and increased hydrothermal activity in the crust sampled by receiver functions recorded at BACH stations would also elevate v p/v s (O’Connell and Budiansky, 1974). With only these data we can not definitely resolve the presence of partial melt; however, we surmise that it is not present in large amounts. Finally, recent InSAR (interferometric synthetic aperture radar) observations in the Changbaishan region indicate a broad region of rapid subsidence (possibly up to 9 cm/year) essentially coincident with the lateral extent of the low velocity anomaly inferred in this study (Kim et al., 2001). This may imply that the hot region in the midcrust is
cooling rapidly due to circulation of meteoric water in a broad region of heavily cracked rock.
Acknowledgments This was a cooperative project between the Research Center of Geophysical Exploration, Chinese Seismological Bureau, and the State University of New York at Binghamton, covered under the Sino-US Seismological Research Protocol signed between USGS/NSF and CSB in 1980. We would like to thank IRIS/PASSCAL for providing the seismometers, as well as software, technical assistance and field support. Installation of the broadband network would not have been possible without the field assistance of Noel Barstow of PASSCAL and IRIS intern Margaret Boettcher. PASSCAL software assistance was provided by Sid Hellman and Paul Friberg. We thank Charles Ammon and Malcolm Sambridge for providing software used in this study. Most figures and maps were produced using GMT (Wessel and Smith, 1991). Careful reviews by F. Darbyshire, L. Pujades and editor K. Furlong are greatly appreciated, and EAH thanks R. van der Hilst and B. Hager for assistance and advice. The coordination and cooperation of the Changchun Office of the CSB and the staff of the RCGE, under the leadership of Dr. Zhang, was essential for the installation, operation and success of the temporary network. This project was funded under DOE grant DE-FC04-98A178918.
References Ammon, C.J., 1991. The isolation of receiver effects from teleseismic P waveforms. Bull. Seismol. Soc. Am. 81, 2504 – 2510. Ammon, C.J., Randall, G.E., Zandt, G., 1990. On the nonuniqueness of receiver function inversions. J. Geophys. Res. 95, 15303 – 15318. Brantley, S.L., Evans, B., Hickman, S.H., Crerar, D.A., 1990. Healing of microcracks in quartz: implications for fluid flow. Geology 18, 136 – 139. Chevrot, S., van der Hilst, R.D., 2000. The Poisson ratio of the Australian crust: geological and geophysical implications. Earth Planet. Sci. Lett. 183, 121 – 132. Christensen, N.I., 1996. Poisson’s ratio and crustal seismology. J. Geophys. Res. 101, 3139 – 3156. Christensen, N.I., Fountain, D.M., 1975. Constitution of the lower continental crust based on experimental studies of
E.A. Hetland et al. / Tectonophysics 386 (2004) 157–175 seismic velocities in granulite. Geol. Soc. Amer. Bull. 86, 227 – 236. Cui, Z., Wei, H., Liu, R., 1995. A textural research of record history about eruption activity of the Tianchi volcano. Volcanic Active and the Human Environment. Seismological Press, Beijing, pp. 36 – 39 (in Chinese). Darbyshire, F.A., Priestley, K.F., White, R.S., Stefa´nsson, R., Gudmundsson, G.B., Jakobsdo´ttir, S.S., 2000. Crustal structure of central and northern Iceland from analysis of teleseismic receiver functions. Geophys. J. Int. 143, 163 – 184. Du, Z.J., Foulger, G.R., 2001. Variations in crustal structure across central Iceland. Geophys. J. Int. 145, 246 – 264. Fan, Q., Liu, R., Li, D., Li, Q., 1999. Significance of K–Ar age of bimodal volcanic rocks at Wangtian volcano, Changbaishan Area. Chin. Sci. Bull. 44, 660 – 664 (in Chinese). Fountain, D.M., Christensen, N.I., 1989. Composition of the continental crust and upper mantle; a review. Geol. Soc. Amer. Bull. 172, 711 – 742. Guo, L.C., Ma, S.Z., Zhang, Y.S., 1996. The magma study of Changbai volcano by CT technique. CT Theory Appl. 5, 47 – 52 (in Chinese). Jin, B., Wei, X., 1994. Researching Volcanic Geology in Changbai Mountain. Northeast Korean Nationality Education Press, Yanbian (in Chinese). Julia`, J., Vila, J., Macia`, R., 1998. The receiver structure beneath the Ebro Basin, Iberian Peninsula. Bull. Seismol. Soc. Am. 88, 1538 – 1547. 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. Kennett, B.L.N., Engdahl, E.R., 1991. Travel times for global earthquake location and phase identification. Geophys. J. Int. 105, 429 – 465. Kim, S.W., Won, J.S., Kim, J.W., Moon, W.M., Min, K.D., 2001. Multi temporal JERS-1 SAR investigation of Mt. Baekdu stratovolcano using differential interferometry. Geosci. J. (Seoul) 5, 301 – 312. Langston, C.A., 1979. Structure under Mount Rainier, Washington, inferred from teleseismic body waves. J. Geophys. Res. 84, 4749 – 4762. Liu, J., 1983. Study on Cenozoic Volcanic Activity in Changbaishan Area. Beijing Science and Technology Press, pp. 343 (in Chinese). Liu, R., Chiu, S.H., Cai, L.Z., Wei, H.Q., Yang, Q.F., Xi, Z.Q., Fu, G.C., Zhong, J., 1997. Chronological study and its implications of recent major eruptions of Tianchi volcano in Changbai Mountain. Sci. China 27, 437 – 441 (in Chinese). Liu, R.X., Wei, H.Q., Lee, J.T., 1998. The Recent Eruptions of the Tianchi Volcano in the Changbaishan. Science Publishers, Beijing, China. (in Chinese). Mavko, G.M., 1980. Velocity and attenuation in partially molten rocks. Geophys. J. Int. 85, 5173 – 5189. O’Connell, R.J., Budiansky, B., 1974. Seismic velocities in dry and saturated cracked solids. J. Geophys. Res. 79, 5412 – 5426.
175
Owens, T.J., Zandt, G., 1997. Implications of crustal property variations for models of Tibetan plateau evolution. Nature 387, 37 – 43. Phinney, R.A., 1964. Structure of the Earth’s crust from spectral behaviour of long-period body waves. J. Geophys. Res. 69, 2997 – 3017. Sambridge, M., 1999. Geophysical inversion with a neighborhood algorithm: i. Searching a parameter space. Geophys. J. Int. 138, 479 – 494. Shangguan, Z., Sun, M., 1997. Mantle-derived rare-gas releasing features at the Tianchi volcanic area, Changbaishan mountains. Chin. Sci. Bull. 42, 768 – 771. Shangguan, Z., Zheng, Y., Dong, J., 1997. Material sources of escaped gases from the Tianchi volcanic geothermal area, Changbai Mountains. Sci. China 40, 390 – 397. Sheriff, R.E., Geldart, L.P., 1982. Exploration Seismology: Volume 1. History, Theory and Data Acquisition. Cambridge University Press, Cambridge. Shin, J.S., Baag, C.E., 2000. Moho depths in the border region between northern Korea and northeastern China from waveform analysis of teleseismic pMP and pP phases. Geosci. J. (Seoul) 4, 313 – 320. Song, J.L., Wu, F.T., Hetland, E.A., Zhang, X.K., 2003. Wide angle refraction/reflection P wave velocity structure of Changbai Volcano, China. In preparation. Tang, J., 1998. The electronic structure of Tianchi volcano area by magnetotelluric. In: Liu, R.X., Wei, H.Q., Lee, J.T. (Eds.), The Recent Eruptions of the Tianchi Volcano in the Changbaishan. Science Publishers, Beijing, China, pp. 108 – 123 (in Chinese). van der Hilst, R.D., Kennett, B.L.N., Shibutani, T., 1998. Upper mantle structure beneath Australia from portable array deployments. In: Braun, J. (Ed.), Structure and Evolution of the Australian Continent, volume 26 of Geodynamics. American Geophysical Union, Washington, DC, pp. 39 – 57. Wessel, P., Smith, W.H.F., 1991. Free software helps map and display data. EOS Transactions, American Geophysical Union 72, 441, 445–446. Wilkens, R., Simmons, G., Caruso, L., 1984. The ratio of vp/vs as a discriminant of composition for siliceous limestones. Geophysics 49, 1850 – 1860. Zhang, J., Langston, C.A., 1995. Dipping structure under Doubes, Belgium, determined by receiver function modeling. Bull. Seismol. Soc. Am. 85, 254 – 268. Zhang, C.K., Zhang, X.K., Zhao, J.R., Liu, B.F., Zhang, J.S., Yang, Z.X., Hai, Y., Sun, G.W., 2002. Crust–mantle structure of the Changbaishan Tianchi volcanic region and its vicinity: an exploratory study and inferences. Chin. J. Geophys. 45, 862 – 871. Zhu, L., Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver functions. J. Geophys. Res. 105, 2969 – 2980.