Crust Macrofracturing as the Evidence of the Last Deglaciation

Machine learning methods were applied to reconsider the results of several passive seismic experiments in Finland. We created datasets from different stages of the receiver function technique and processed them with one of basic machine learning algorithms. All the results were obtained uniformly with the $k$-nearest neighbors algorithm. The first result is the Moho depth map of the region. Another result is the delineation of the near-surface low $S$-wave velocity layer. There are three such areas in the Northern, Southern, and central parts of the region. The low $S$-wave velocity in the Northern and Southern areas can be linked to the geological structure. However, we attribute the central low $S$-wave velocity area to a large number of water-saturated cracks in the upper 1-5 km. Analysis of the structure of this area leads us to the conclusion that macrofracturing was caused by the last deglaciation.


Introduction
Passive seismic experiments SVEKALAPKO (Bock et al., 2001) and POLENET/LAPNET (Kozlovskaya et al., 2006) spurred a large number of research activities that targeted the lithosphere structure of the northern part of the Baltic Shield (Bruneton et al., 2002;Alinaghi et al., 2003;Aleshin et al., 2006;Vecsey et al., 2007;Kozlovskaya et al., 2008;Uski et al., 2011;Grad and Tiira, 2012;Silvennoinen et al., 2014;Vinnik et al., 2016).A significant part of these studies is based on the receiver function technique (Alinaghi et al., 2003;Aleshin et al., 2006;Kozlovskaya et al., 2008;Grad and Tiira, 2012;Vinnik et al., 2016).With data from a wide network of seismic stations across the region, this method yielded several important results on the structure of the mantle of northern and southern Finland (Alinaghi et al., 2003;Frassetto and Thybo, 2013;Vinnik et al., 2016), and played a key role in the implementation of a three-dimensional seismic model of central and southern Finland (Kozlovskaya et al., 2008).An important structural feature of the region is the presence of a thin surface layer with low shear wave velocities Vs.The layer thickness varies around 1 km, the Vs values are down to 40 percent less than the average over the crust (Aleshin et al., 2006).The presence of this layer was first noted in (Pedersen and Campillo, 1991;Grad and Luosto, 1992;Grad and LUOSTO, 1994) in the study of the seismic surface wave attenuation.They reveal not only low values of -wave velocity in ORCID(s): 0000-0001-7489-2951 ( I.M. Aleshin); 0000-0003-0324-9795 ( K.I.Kholodkov) the shallow crust but also low values of seismic quality factor (Q-factor), which increase rapidly down to the depth of about 1 km.According to Grad et al. (1998), the Q-factor is suitable to differentiate the Archean and the Proterozoic crust of the Fennoscandian Shield.In Grad and Luosto (1992) showed that the increased density of cracks provide for lower Q-factor for the top 1 km of the crust.This explains the low -wave velocity.The presence of the low -wave velocity surface layer (LVSL) was confirmed later in Aleshin et al. (2006).However, Kozlovskaya et al. (2008) have shown that the presence of such a layer is independent of the age of the rocks, it is present in a significant part of both Archean and Proterozoic regions.The research of the Moho shape adopts the results of our previous studies in the SVEKALAPKO and POLENET/LAPNET projects (Aleshin et al., 2006;Kozlovskaya et al., 2008;Silvennoinen et al., 2014) supplemented by data from earlier research (Dricker et al., 1996;Aleshin et al., 2019) (see Code Availability section for this data).Additionally, we adopted 1D profiles to construct the 3D seismic image in the arias of interest.
We had a relatively small amount of input data.The spline interpolation like in (Horspool et al., 2006) or kriging technique as described in (Kozlovskaya et al., 2008) can not be considered optimal in this kind of problem.Instead, we preferred to use a machine-learning approach and chose the k-nearest neighbors (kNN) algorithm for its simplicity.
One more reason to choose kNN is the data behind the analysis as it consists of separate sets of three different data types.The first one is the "yes-or-no" presence of the converted phase right after the -wave arrival in the receiver function's waveform.The second type is a series of -wave velocity vs. depth pairs.The third type is crust thickness beneath the stations.Universally, the machine-learning approach allows us to process all of these types.In the "Data" and "Method" sections we described the input dataset and covered the processing algorithm features, correspondingly.
We applied the method to the Moho depth data in the section "Moho Boundary Shape".Notably, the results appeared to be in good alignment with previous works on the Moho shape of this particular region.In the next section, we applied the machine learning approach to locate and describe the structure of the LVSL regions.In the "Discussion" section we argue in favor of post-glacial rebound as the origin of the low -velocity layer in central Finland.

Data
For our analysis, we predominantly utilized the results of the receiver function study of passive seismic experiments data, namely SVEKALAPKO and POLENET/LAPNET.Additionally, we used the results of the previous receiver function study for stations LVZ and APA (Dricker et al., 1996); PITK and KEMI (Aleshin et al., 2019).Finally, we turned to account for the models derived from wide-angle seismic experiments (Janik et al., 2007(Janik et al., , 2009)).
The LVSL under the station is manifested by the strong converted phase arriving within the first seconds after the main -wave.Such a phase is clearly seen on the receiver function waveform (Figure 1, panel A, the FH03 station).
We used it to construct the first dataset (dataset "A") that represents the spatial distribution of the LVSL.Dataset members consist of the station name, location, and logical condition of whether the LVSL is present or not.A few additional dataset members were included from the models derived from wide-angle seismic experiments (Janik et al., 2007(Janik et al., , 2009)).It was not possible to draw a definite conclusion on the absence or existence of the LVSL under the stations FB11, FD07, LP11, FJ01, LP61, ARE0, LP75.These points were excluded from the dataset.As a result of the analysis, we obtained a set of points with the given values of the Boolean type.The dataset "A" enables us to compute the probability of the LVSL presence at an arbitrary point.To get more details on the structure we require a spatial distribution of elastic parameters.We have obtained it from the 1D -velocity profiles by joint inversion of -receiver functions and surface wave dispersion.Those profiles were adopted from (Kozlovskaya et al., 2008).
For the sake of convenience, we quantized the original layered models with 0.1 km step.The choice of the step is more or less volitional.The lower step limit is determined by computation time and complexity.The upper step boundary is limited to the thickness of the thinnest layer across all stations, which appeared to be 0.7 km for our case.
Therefore, we've got a second dataset, "B".It contains the station name, location (latitude, longitude?and depth), and the corresponding -velocity.Members of the dataset are irregularly spatially distributed.1D -velocity profile also provides for determining the Moho depth under the station (see fig. 1, panel B).POLENET/LAPNET (Silvennoinen et al., 2014) and SVEKALAPKO (Kozlovskaya et al., 2008) projects used the same approach to determine the Moho depth.Additionally, the Moho depth for stations LVZ, APA, and PETK, KEMI were adopted from (Dricker et al., 1996) and (Aleshin et al., 2019), respectively.Numerical values of the Moho depth for all the stations are unified into the third type dataset "C".We used this dataset to construct the Moho depth map.

Method
At this point, we have three separate datasets with different types of data: two-dimensional logical values (dataset "A"), three-dimensional floating-point values (dataset "B"), and two-dimensional floating-point values (dataset "C").
We need to obtain the synthetic value for all types of data in datasets at an arbitrary location and depth, where applicable, we opted for using machine-learning methods to do this universally.To perform this task we have to, firstly, define the distance rules for each type of data.For data in dataset "A" and "C", we can utilize the orthodrome (the geodetic line length).However, due to three-dimensional data in the dataset "B", the cartesian distance is more convenient.
Because the region of interest is quite small we switched from geographical latitude and longitude to a local Cartesian coordinate system ( , ) using approximate formulae = 2 ( − 1 ) cos ∕360; = 2 ( − 1 )∕360; = ( 2 − 1 )∕2; (1) Here is the Earth's radius, 1 , 2 are the minimal and the maximal latitudinal boundaries of the region, correspondingly; 1 is the minimal longitude.Also, the station elevation was set to be zero as we did not include the terrain features into the consideration.There is no prominent influence of the simplifications on the conclusion; we used that due to convenience only.Thus, the distance between two points ⃗ = { , , } and ⃗ = { , , } was defined as -axis is assumed normal to the horizontal plane ( , ) plane and zeroes at sea level.The meaning of the scaling factor will be explained later.
As soon as we have defined the distance we should choose a formal mathematical model (algorithm) to be able to calculate the target value for any given location.In general, the algorithm contains a few inner tuning parameters.We can fix them with the input data (the algorithm learning).We've chosen the kNN-algorithm, one of the basic machinelearning algorithms.The kNN does not imply any preliminary training procedure.It is evaluated only to predict the target value at the given location ⃗ by the set of known values Here ( ) is the weight function that depends on the distance from the current point ⃗ to the location of the -th dataset member.The summation is performed over the dataset members closest to the point ⃗ .Usually, the weight function is chosen in the form of a power law: In this paper, we used this law with the value = 1.The kNN -algorithm is called a lazy learner because it implies no separate training process.In this algorithm, the dataset is stored in the memory in a certain spatial order based on the distance (2).Formula (3) shows that inner parameters of the algorithm are just points of the input data nearest to a given position.The computational consumption originates from the search of the members closest to the given point.
The complexity of the brute force selection from elements is about ( ), which is quite inefficient.However, if the data is stored using the binary space partitioning (BSP) scheme the complexity can be reduced to (log ).In our work, we used NumPy (Harris et al., 2020) arrays to store the data, which enabled hassle-free k-d tree-like BSP or ball-tree storage organization.
Along with inner parameters, the algorithm used in machine-learning often has hyperparameters, which do not change during training.In our case, these hyperparameters are the number of nearest neighbors and scale factor from (2).We chose to determine the optimal values of (and for dataset "B") using the cross-validation technique (Hastie et al., 2009).The main idea is to split the input data into training and validation blocks.Training blocks are used to fit the internal parameters of the model.Then, the trained model is used to predict the data of the validation block.The comparison of the predicted data and the validation data is used to estimate the training quality.
Within the actual implementation, the initial set is to split into blocks, which we will call folds, each of the same size.Each fold is consistently served as a validation block, and the remaining − 1 folds act as the training blocks.The main advantage of this strategy is that each piece of data is used both for training and for quality control.There are a lot of strategies for splitting data.For dataset B we chose = 5 (see below).In problems with a lack of data, that is the case for datasets "A" and "C", the value of is appropriate to be chosen equal to the number of members of the source dataset so that each fold contains just a single point.This strategy is called Leave-on-Out.To estimate the training quality, we have to define the quality function.In regression analysis one of the simplest choices is a mean squared difference between the values calculated from the training block and the corresponding samples in the validation block.
Essentially, this metric is a measure of the mean interpolation error.We have chosen such a quality function to determine hyperparameters or , in the Moho boundary shape (dataset "C") and -velocity 3D-model problems (dataset "B"), correspondingly.However, this metric is unsuitable for the calculation of areas covered with a lowvelocity layer (dataset "A").For this problem, we have to formulate the task as a binary classification problem.Within this classification, we divide all the points into two classes.Let the point belong to class 1 if the low -velocity layer at the location of the point exists and to class 0 otherwise.Therefore, our goal is to compute the probabilities that a given point on the surface belongs to class 1 or class 0. In this setting, receiver operating characteristic (ROC) is commonly employed.The ROC curve is created by plotting the true positive rate (probability of detection) versus the false positive rate (probability of false alarm) generated by the binary classifier.The area under the curve (AUC) gives a measure of the classifier quality (more detail can be read elsewhere, (Fawcett, 2006)).The rest of the calculation details are the dataset specific and will be described later.

Moho Boundary Shape
For this region, the depth to the seismic Mohorovicic discontinuity approximately equals the compositional crust- Several reasons can explain this difference.In the controlled source approach the Moho depth is obtained by -waves in most cases, while the receiver function utilizes -or -waves from teleseismic sources converted at the Moho.Also, the frequency content of the signal from teleseismic source data is rather distinct from that in near-vertical reflection, wide-angle reflection, and refraction experiments.Furthermore, rays in both the near-vertical setup and the receiver function are sub-vertical apart from the wide-angle approach.The other reasons are region-specific.In the southern part of the region, the lowermost part of the crust has a high-velocity layer (Janik et al., 2007) that masks the high-velocity contrast at the Moho for -wave (Janik et al., 2009).Also, there is a tangible seismic anisotropy in the upper mantle (Vecsey et al., 2007;Vinnik et al., 2016;Yanovskaya et al., 2019) that causes the discrepancy in sub-lateral and sub-vertical -wave travel times.The number of nearest neighbors is only one hyperparameter in the Moho boundary shape problem.To determine an optimal value of we adopted the Leave-one-Out cross-validation procedure.We used the quality function (5) (we set = 0).The optimal value of = 4 minimizes the functional (see fig. 2).The corresponding mean interpolation error is 3.7 km.The map in Fig. 3 depicts the depth to the Moho boundary, red --data from this study.The map was constructed with formulas (1),( 2),(3) where = 4.The maximum thickness of the crust appears to be in the center of the southern part of the region which corresponds to the results described in Alinaghi et al. (2003); Kozlovskaya et al. (2008); Silvennoinen et al. (2014).Additionally, there is a correlation of the Moho boundary with the Central Finnish Granitoid Complex shapes, which is also discussed in cited papers.Fig. 4 shows the -profile along FD13-FJ01.The black line marks the Moho depth obtained in Kozlovskaya et al. (2008) with the kriging technique.The solid red line denotes the Moho boundary obtained in this study with the kNN approach.Both lines agree well, showing good consistency among the approaches.Dashed red lines show interpolation error intervals.The map (Fig. 3) and the profile (Fig. 4) reveal the crust thickness anomaly under station FG01.The same was noted by Kozlovskaya et al. (2008).It is worth noting that interpolation misfit surges in the vicinity of FG01, hitting values times the mean error.Thus, we lack the confidence to conclude on the Moho depth there.

Map of the low -velocity Layer in the uppermost crust
Now, we can apply the machine learning approach to dataset "A" to build a map of the LVSL coverage in the region.
At first, we have to determine the optimal number of nearest neighbors .We used the ROC-AUC metric and applied the Leave-one-Out strategy for cross-validation, similar to the Moho shape calculation.The larger ROC-AUC value corresponds to the better classification model.From the left panel of 5, one can see that the optimal number of nearest neighbors is reached at = 4.The shape of ROC at this value is plotted on the right panel.To construct the map of LVSL from the dataset "A" we substituted = 0 and = 4 in formulae 1, 2. The calculation was performed at nodes ⃗ of the regular surface grid.By the definition, every member of the input dataset "A" is either 0 or 1.It leads to the function = ( , ⃗ ) from formula ( 5) is constrained to 0 ≤ ≤ 1.Thus might be treated as the probability of the presence of the LVSL at the given point ⃗ .To perform the final classification we had chosen two threshold values 0 and 1 .Whether the was less than 0 = 0.4 we concluded that there was no LVSL at the point ⃗ .As a result, the point was classified as class "0".Otherwise, when exceeded 1 = 0.6 we decided that there was an LVSL at this pint, and ⃗ was classified as class "1".The interim values of that comply with ( 0 < < 1 ) remained unclassified and formed the "buffer zones" between two classes.
The map obtained is shown in Fig. 6.The map contains three separate LVSL areas in the northern, central, and southern parts of Finland.Stations excluded from the analysis are mostly located near the borders of those areas like FB11, FD07, LP61, or near the outer border of the region like ARE0, LP75.Such distribution explains partly the difficulties of the initial classification.With new results, we can suspect the existence of the LVSL under the station LP61.There appears to be no such layer under the stations LP11, FJ01, and in point HT.The stations FD07 and FB11 are located inside the "buffer zones" so we can not decide the LVSL under these points.The southern and northern LVSL areas could be linked to the geology: the metamorphosis of bedrock and the presence of Paleozoic sedimentary rocks atop the archean basement, respectively.In the central part of Finland, the LSLV can be explained with neither the composition nor age of the rocks.As seen in Fig. 6 the LVSL area in central Finland covers both the Archean and Proterozoic rocks.The origin of the layer should be discussed separately.But before that, we used dataset "B" to highlight the structure under this area.We have to mention that when we quantized the 1D -velocity profiles we got vertical distances between the points of the same profile much smaller than distances between the stations.Here, to illustrate this issue we would stick with the alternative formulation of kNN where the term "nearest neighbors" assumes all the points inside the sphere of radius R around the given position, see Fig. 7, left.This corresponds to set = 1 in the formula (2).With moderate R the nearest points would originate from the same profile.The issue could be solved by using the ellipsoid instead of the sphere and setting > 1.We used a fixed the nearest neighbor number instead of the radius of the sphere and the factor controls the -axis scaling.To determine the optimal and we used the cross-validation technique for dataset "C".The dataset "C" has quite a large number of points so we've split it into 5 folds of approximately 8000 points each.The quality function ( , ) calculated for = 1, 2, ..., 150 and = 1, 2, ..., 150 is shown on Fig. 7, right.Such a distribution is intrinsic to concurrent optimization problems, as seen in, for example, Chiu et al. (2009).Optimal pairs of and are located where the function ( , ) is minimal.
These values of are placed at the axes origin on Fig. 7, right.However, this isn't the case in our problem because picking such low values will lead to very low K and and subsequently to an overinterpretation.Larger values of and correspond to smoother interpolation results, and to larger ( , ).Meanwhile, the error of interpretation cannot exceed the error of the receiver function method, which is typically about 0.1 km/s.Overcauteously, we chose 0.05 km/s for the ( , ).This level provides for the optimal pairs of and .For example, picking points A, B, and C from Fig. 7, right net visually equivalent results.So we picked point B, which is the closest to the origin.Such a choice is obvious for multi-objective optimization problems Chiu et al. (2009).We visualized the horizontal section of the -wave velocity at the surface level ( = 0), see Fig. 8, left.The low velocity region matches the probabilities computed with dataset "A" well.There are strong -velocity variations of the LVSL in archean, proterozoic, and transition regions.We can speculate that the thickness of the LVSL could depend on rock composition, however, investigation of this dependency is out of the scope of the current work.We also imply that the -velocity depends on the amount of rock fracturing and water saturation.

Discussion
In our study, we do not make a particular interpretation of -wave velocity values of the low -velocity layer in the uppermost crust.The most interesting fact is that this layer has been revealed both in Archean and Proterozoic areas and it is existing in bedrock units with different compositions.Moreover, it is not correlating with the rock units composed of low-velocity rocks, like granites.For example, the layer is not detected in the Central Finland Granitoid Complex.
Therefore, its origin seems to be not dependent on the tectonic age or composition of the bedrock and rock fracturing is the main factor responsible for low velocities.On the other hand, this layer is not distributed uniformly around the Fennoscandia and it is present in several large-scale areas indicated by dark blue in Fig. 8.These two features (e.g.absence of correlation with tectonic age and rock composition and systematic distribution in large-scale areas) need to be explained.In our study, we used the converted waves with dominating periods of 2 s, which means that the wave is averaging fractures of different scales up to 1 km.Such large-scale fractures reaching the depth of 1 km could not form due to rock weathering, but rather due to some large-scale processes.One of the explanations of this layer could be that it has been formed due to the processes related to multiple episodes of deglaciation of the Fennoscandian Shield.
Relation between intraplate seismicity and processes of deglaciation (Post-Glacial Rebound, PGR) has been studied previously by numerous authors, see Turcotte and Schubert (2002); Lagerbäck and Sundh (2008); Slunga (1991).The seismicity at the ice sheet margins is attributed to unloading processes and to ice melting (Wu and Johnston (2000)).
According to Spada et al. (1991), the PGR is the valid reasoning for normal and lateral crustal deformation at polar regions.The PGR appeared to have a solid role in stress redistribution.[Steffen et al., 2012;Steffen et al., 2014] unfold that in the cases of Greenland and Antarctica the earthquakes mainly credit to the PGR.However, the previous studies of seismicity in the Fennoscandia (Lagerbäck and Sundh, 2008;Slunga, 1991;Uski et al., 2003Uski et al., , 2012) ) were mainly trying to explain the origin of moderate-to-large magnitude earthquakes (M>3) with sources corresponding to large-scale faults.Studying the effect of the PGR on the formation of smaller scale fractures in the bedrock has not been done previously.
Nowadays climate change caused rapid warming and melting of ice sheets of Greenland and Antarctica.Before this massive melting, these areas were considered as generally aseismic, which has been explained by the suppressing of earthquakes due to ice sheet load.The rapid deglaciation results in significant stress change in these areas, both in the crust and in the mantle, and hence in new seismic events in these areas.That is why more detailed studies of the connection between PGR and bedrock fracturing at different scales are possible nowadays.In particular, Olivieri and Spada (2015) studied the relationship between present-day crustal deformations and spatio-temporal distribution of recent earthquakes in Greenland.They showed that nowadays, Greenland undergoes the significant deformation due to two processes: the melting of the Greenland Ice Sheet (the "elastic rebound") and the PGR from late-Pleistocene.So, Olivieri and Spada (2015) suggested two distinct mechanisms responsible for seismicity: regional crustal earthquakes are attributed to the PGR, and local events are attributed to the ER and present-day ice melting.In the areas with the largest ice sheet thickness, the seismicity is suppressed nowadays.If we consider deglaciation of Greenland as a modern analog of processes of Fennoscandia deglaciation, we may speculate that the LVSL in the Fennoscandia was formed due to ER in the areas which underwent rapid deglaciation in the past.Previous studies of Fennoscandian ice sheet melting and retreat models and reconstructions have been summarized by Stroeven et al. (2016).They showed that the deglaciation rate was significantly different across the area investigated in our study.These studies are based on geomorphological and geochronological data.However, the amount of these data is limited, which increases the uncertainty in the modeling of ice sheet deglaciation processes.We think that our study provides new, previously not considered evidence of the effect of deglaciation on the fracturing of the uppermost bedrock that could be used both for upgrading of Fennoscandia deglaciation models and in predicting the consequences of climate change and melting of Greenland and Antarctica Ice Sheets.

Conclusions
The results obtained continue and generalize the series of studies conducted by the authors.This time, we performed a joint interpretation of the results of previous studies supplemented by data for several new stations using a modern mathematical apparatus.The results showed the effectiveness of universal machine learning methods for analysis and generalization.Some algorithms of machine learning benefit highly from a data deficit that is typical in many geophysical studies.Further progress in exploring the region is hindered by the lack of data on a significant part of the Russian territory.

Acknowledgments
The majority of this work would not be possible without our late colleague Dr. Grigory Kosarev.

List of Figures
mantle boundary.The shape of the Moho boundary with receiver function technique was previously studied for several areas in the region: Kola inKosarev et al. (1987), Southern Finland inAlinaghi et al. (2003);Kozlovskaya et al. (2008);Grad and Tiira (2012), Northern Finland inSilvennoinen et al. (2014).In a series of works the Moho was investigated with controlled source seismic experimentsGrad and Luosto (1992);Janik et al. (2007);Uski et al. (2011).InGrad and Tiira (2012);Uski et al. (2011) the crust thickness was studied as part of the European lithosphere structure.In the paper bySilvennoinen et al. (2014), the topography of the Moho under the northern part of Finland was drawn by a joint interpretation of the receiver function and controlled source experiments.In our paper, however, we only use the data obtained with the receiver function technique.We opted for this limitation because as shown inSilvennoinen et al. (2014) the Moho boundary estimated from the receiver function may differ from that of a controlled source experiment.

Figure 1 :
Figure 1: A. Receiver function waveforms of two SVEKALAPKO stations, FH03 (thin line) and FA05 (thick line).FH03 waveform shows clear arrival of the converted phase near the first second after the -wave arrival.This indicates the -velocity jump immediately beneath the station.FA05 doesn't show such a phase.The letters M 1 and M 2 pick to the Moho border phase.B. Corresponding 1D -velocity profiles obtained through the receiver function inversion(Kozlovskaya et al., 2008).The Moho borders are pointed with the letters M 1 and M 2 .

Figure 2 :
Figure 2: The plot of the quality function (5) versus the number of nearest neighbors derived from the Moho boundary depth dataset.

Figure 4 :
Figure 4: Depth to Moho along with the FD13-FJ01 profile.Black and red lines mark a cross-section of the Moho shape obtained with the kriging technique and the machine learning kNN method, correspondingly; red dashed lines show the kNN interpolation errors.The line oscillations at distances of 550-600 km are likely due to erratic data for station FG01.

Figure 5 :
Figure 5: The plot of the ROC-AUC versus the number of nearest neighbors (left).The plot of ROC for the optimal value of = 4 (right).

Figure 6 :
Figure 6: Map of the low-velocity surface layer calculated from the dataset "A".Areas in which the probability of the presence of the LVSL exceeds 60% are colored in blue.Olive shaded areas denote the probability is less than 40%.Triangles indicate the position of seismic stations and squares indicate the points from wide-angle experiments.They are white or black, depending on the presence or absence of the LVSL at the point, respectively.Initially excluded points are marked red.

Figure 7 :
Figure 7: Left: Illustration of the role of scale factor.Right: plot of root-mean-square deviation (formula (5)) for the different values of and .

Figure 8 :
Figure 8: Left: Velocity map of -waves at the surface.Plus-marked ("+") area matches the central LVSL-area on Fig. 6.Right: vertical section of the -wave velocity of AB profile.