Crustal velocity structure under the RUKSA seismic array (Karelia, Russia)

Construction of a Velocity Structure by the Simulated Annealing Method

[13] We modeled the medium under the array by a set of horizontal layers underlain by a homogeneous half-space. The model parameters m included P and S wave velocities, density, and layer thicknesses. The model was constructed through the minimization of the functional

+ w₃ |t^obs₄₁₀ - t^syn₄₁₀|.

A synthetic receiver function Q_syn( m, t) was calculated by the formula [Kind et al., 1995]

The response components of the layered medium H_Q ( m, w) and H_L ( m, w) were calculated by the Thomson-Haskell method [Haskell, 1962]. ( R(w) is the surface wave dispersion curve, and w_a are empirically chosen weights.)

[14] The minimization problem formulated above is strongly nonlinear and ill-conditioned. The traditional method used for solving problems of this type is based on the regularization method [Tikhonov and Arsenin, 1979]. The minimization reduces to the solution of a system of linear equations that can be found very rapidly. However, the model obtained from an implementation of this method [Kosarev et al., 1987] is essentially dependent on the initial approximation whose choice is generally complicated.

[15] Statistical (nongradient) methods (primarily, the Monte Carlo method) "scan" the entire space of parameters in the process of solution and are free from the necessity of specifying the initial approximation. However, the direct application of statistical methods to the search for the solution in a multidimensional space involves time consuming computations, which makes the practical use of these methods unrealistic. A possible way of solving this problem is to seek the solution in three stages: identification of one or more promising regions in the initial space of potential models, careful examination of each of these regions, and "fine adjustment" of the best solution inferred at the preceding stage. Such an approach is realized in the modern modification of the simulated annealing (SA) method [Ingber, 1989]. The software implementation of this method reduces the time of the search for the solution to an acceptable value (no more than one hour on a standard PC for the present case).

[16] To invert data, we used a model consisting of 11 layers on a half-space. Velocities of S waves in the layers and the half-space and thicknesses of some layers were varied, whereas the ratio V_p/V_s was fixed (1.73 in the crust and 1.8 in the mantle) and the density was determined by the Birch formula. The chosen time interval of the receiver function minimization ( - 2 s to 8 s) bounds the overall depth of the model by about 70-80 km. Deeper layers have no effect on the receiver function. However, in order to use traveltimes of waves converted at the 410-km boundary, the velocity structure should be specified to a depth of 410 km. For this purpose, our current 11-layer model (more specifically, the line defining the velocity in the half-space) was continued downward until its intersection with the IASP91 standard model of the Earth. IASP91 velocities were used for the calculation of tps410 in the depth interval from the intersection point to 410 km.

Figure 4

[17] The SA method does not require the specification of an initial approximation. However, it is necessary to define intervals within which the variable parameters can be changed. The resulting model must lie within these intervals. The variation limits of P and S wave velocities and the inversion results are shown in Figure 4. The P velocity variation limits were chosen in such a way that they included the whole range of velocity values published in the recent literature for the Baltic Shield. For example, the P velocity in the layer at depths of 10-15 km (the first layer of the middle crust) was varied from 6.0 km s ^-1 to 6.4 km s ^-1 in agreement with the seismic velocity structures obtained for this region [Goncharov et al., 1991]. The upper part of the model consists of six layers. The thickness of the two uppermost layers was varied from 0.1 km to 1 km, whereas the thickness of the four remaining layers was fixed (2 km). Such a detailed specification of the upper part of the model was due to the necessity of fitting the intense phase of the first second. The remaining part of the crust includes five layers of variable thicknesses that could change from 5 km to 10 km. the mantle is described by a half-space. During optimization, 6000 to 20,000 models were generated.

Figure 5

[18] To estimate the effect of each type of the observed data on the resulting model, we performed the following experiment (Figure 5). The receiver function alone was first inverted. As a result, we obtained a model shown by the red line in Figure 5. Then the receiver function was supplemented by either surface wave dispersion data (blue line) or the traveltime of converted waves from the 410-km boundary (green line). Finally, the fourth curve (black line) represents a model obtained by using the whole data set. Comparison of models reveals that the data of the receiver function alone underestimate the velocities in both the crust and the mantle. The receiver function data supplemented by either the data of surface waves alone or by the tps410 traveltime data alone yield virtually the same and, apparently, correct velocities in the middle crust at depths of 8 km to 20 km. The surface wave data underestimate velocities in the lower crust and overestimate the mantle velocities. On the contrary, the tps410 data overestimate the lower crust velocities but yield "correct" velocities in the mantle. By correct velocities, we tentatively mean here the velocities obtained from the set of all data (the black line).