A MATLAB Package for Calculating Partial Derivatives of Surface-Wave Dispersion Curves by a Reduced Delta Matrix Method

Abstract

Various surface-wave exploration methods have become increasingly important tools in investigating the properties of subsurface structures. Inversion of the experimental dispersion curves is generally an indispensable component of these methods. Accurate and reliable calculation of partial derivatives of surface-wave dispersion curves with respect to parameters of subsurface layers is critical to the success of these approaches if the linearized inversion strategies are adopted. Here we present an open-source MATLAB package, named SWPD (Surface Wave Partial Derivative), for modeling surface-wave (both Rayleigh- and Love-wave) dispersion curves (both phase and group velocity) and particularly for computing their partial derivatives with high precision. The package is able to compute partial derivatives of phase velocity and of Love-wave group velocity analytically based on the combined use of the reduced delta matrix theory and the implicit function theorem. For partial derivatives of Rayleigh-wave group velocity, a hemi-analytical method is presented, which analytically calculates all the first-order partial differentiations and approximates the mixed second-order partial differentiation term with a central difference scheme. We provide examples to demonstrate the effectiveness of this package, and demo scripts are also provided for users to reproduce all results of this paper and thus to become familiar with the package as quickly as possible.

1. Introduction

Over the past three decades, a variety of surface-wave exploration methods has become an increasingly important means for inferring the properties of subsurface structures, especially in the near-surface region [1,2]. Compared with the commonly employed seismic body waves, the seismic surface waves (SWs), which primarily include Rayleigh and Love waves in land acquisition, usually propagate in the vicinity of the free surface and exhibit conspicuous dispersion characteristics in seismic records. This dispersion characteristic of SWs is generally delineated by the so-called phase- and/or group-velocity dispersion curves (DCs). The near-surface shear-wave velocities (Vs) can be obtained by inverting these DCs. The inverse problem of DCs, inherently nonlinear, can be solved by using global optimization methods, for instance, genetic algorithm [3], simulated annealing [4], particle swarm optimization [5], etc., or by using various linearized inversion strategies [1,6,7]. The success of either DCs inversion strategy is based on the accurate and efficient computation of surface-wave DCs. In addition, if the latter strategies are chosen, reliable calculations of partial derivatives (PDs) of DCs with respect to the parameters of subsurface strata are indispensable to create the Jacobian matrix (or the sensitivity matrix) in this class of strategies, especially due to the demand of their repeated calculations during the entire inversion process.

The efficient computation of surface-wave DCs benefits not only from the pioneering groundwork of Thomson and Haskell [8,9], but also from the rapid development of computer techniques after the 1960s. Inspired by Thomson and Haskell, numerous approaches have been proposed since then [10,11,12,13,14,15,16,17,18,19,20,21,22,23], particularly to remedy the issue of loss of numerical precision at high frequencies, which is an inherent shortcoming of Haskell’s theory (see reference [23] for more details). Among all these approaches, the one put forward by Dunkin [13] and its improved variants [14,23] are frequently utilized as the foundation of various surface-wave DCs inversion strategies due to their simplicity, efficiency, and effectiveness. Dunkin’s approach is based on the so-called delta matrix theory [13], with which the shortcoming of Haskell’s theory are thoroughly overcome.

Press et al. [10] developed the first computer program for computation of Rayleigh- and Love-wave DCs. In spite of the high efficiency in computing phase-velocity DCs, they applied the numerical differentiation of phase velocity to acquire the group velocity, which is a brute-force, low-precision, and unreliable method. The first category of analytical approaches for calculating the surface-wave group-velocity DCs and their PDs is based on the variational theory [24,25,26,27,28,29]. By virtue of Rayleigh’s principle, Jeffreys [24] first derived expressions for computing the group velocities of surface waves analytically; nonetheless, the PDs of DCs were approximated to the first order accuracy with numerical differentiation. Thereafter, the same idea was adapted to tackle more sophisticated cases, such as Takeuchi et al. [25,26], Anderson [27], Harkrider [28], and Aki et al. [29], just to name a few. On the other hand, to bypass the theoretical complexity of these approaches, Novotný introduced the implicit function method, another analytical approach, first for Love waves [30,31] and then for Rayleigh waves in the case of a layer overlying a half-space [32]. Then Urban et al. [33] succeeded in generalizing this approach to the case of multiple layers for Rayleigh waves. With the help of Knopoff’s algorithm [11,15,16], in combination with the implicit function method, they successfully acquired the analytical expressions for computing the PDs of Rayleigh-wave phase velocities, and the group velocities are also acquired by further using Rodi’s method [34].

To date, the most credit on this subject is given to researchers who focused on the scale of crustal seismology and who were mainly active before the 1990s. However, with the flourishment of the surface-wave methods for near surface characterization [1,2,3,4,5,6,7], this subject has regained new interest in recent years. For example, Cercato recently applied the reduced delta matrix method along with the implicit function theory to achieve the analytical computation of PDs of Rayleigh-wave ellipticity [35] as well as the necessary phase velocity [36]. In order to increase the computational efficiency, he formulated the PDs of Rayleigh-wave phase velocity in a vector-transferring manner instead of the matrix multiplication. Nevertheless, this computational manner makes it intractable to further calculate the PDs of group velocity analytically since more complex second-order partial differentiations must be dealt with. In addition, he only considered the Rayleigh waves in the study. Due to the complexity of this analytical method, it is not widely used in current surface-wave linearized inversion strategies and many researchers have resorted to the global optimization methods as previously mentioned, bypassing the calculation of PDs. However, these methods require a large amount of calculation, particularly when thousands of DCs need to be inverted, as in the context of oil seismic exploration.

Based on the above-mentioned reasons, here we rederive a set of new analytical expressions in the manner of matrix multiplication and present the relevant algorithms to compute the PDs of surface-wave DCs. This work can be taken as another option or a complement to the previous studies. We derive the relevant formulas and algorithms in such another way that is more convenient and more straightforward for considering all of the cases (Love and Rayleigh waves; phase and group velocity). Furthermore, we develop an open-source MATLAB package called SWPD (Surface Waves Partial Derivatives) accompanying this paper, and it has the following merits: (i) it is able to deal with both Rayleigh and Love waves; (ii) it can calculate the PDs of surface-wave phase-velocity DCs analytically or accomplish the calculations for Rayleigh-wave group-velocity DCs by a hemi-analytical method with high precision; (iii) all the relevant expressions and algorithms are derived from a widely accepted frame (i.e., the reduced delta matrix method), which makes the package easy to use and to be extended. It is worth mentioning that the hemi-analytical method is a compromise: it calculates all the relevant first-order PDs analytically and approximates the only second-order mixed PD term with a central finite-difference scheme. By doing so, we still can obtain satisfying results comparable to the analytical approach [32], meanwhile avoiding the intricate analytical derivations for the second-order term.

We first outline the reduced delta matrix theory and its application in modeling the surface-wave phase-velocity DCs. We then present the algorithms and related formulas to calculate analytically the PDs of the phase velocity (Section 3). Next, both the fully analytical and the hemi-analytical approaches for calculating the PDs of the surface-wave group velocity are elaborated (Section 4). After that, we also describe the accompanying MATLAB code with focus on some key scripts or functions (Section 5). Last, four models are utilized for verification of our method and of the package (Section 6).

2. Computations of Surface-Wave Phase-Velocity DCs Based on the Reduced Delta Matrix Method

We first give a brief overview of the reduced delta matrix method [14,23,35] since it constitutes the foundation of the following sections. We strictly follow the formulism of Buchen et al. [23] due to his modern mathematical representations compared with other earlier studies.

This research involved a multilayered earth model (Figure 1), which consists of many homogeneous, isotropic, and elastic layers that overlie a half-space and lie beneath a free surface. Each layer is denoted by the P-wave velocity ${\mathsf{\alpha }}_m$, the S-wave velocity ${\beta }_m$, the density ${\rho }_m$ and the thickness $h_m$. According to the reduced delta matrix method, the surface-wave secular (or dispersion) equation F can be described as an implicit function by the delta matrix recursion:

(1)$\begin{array}{l}{\overline{X}}_{m+1}^{*}={\overline{X}}_m^{*}{\overline{T}}_m^{*}, m=1,2,\cdots ,n\\ F\left(\omega ,c\right)={\overline{X}}_{n+1}^{*}{\overline{V}}_{n+1}^{*}=0\end{array}$

where ω is the angular frequency, c is the phase velocity, m is the sequence number for layers, and n + 1 indicates the half-space. ${\overline{X}}_m^{*}$ is a row vector with different lengths depending on which wave (Love or Rayleigh wave) is concerned. ${\overline{T}}_m^{*}$ and ${\overline{V}}_{n+1}^{*}$ can be regarded as the reduced delta matrices corresponding to the layer m and the half-space, respectively. The elements of these three matrices are listed in Appendix A. In order to obtain the phase-velocity DCs, the frequencies or periods of interest are substituted into (1), followed by a root-finding procedure consisting of the integrated use of the bisection method and the Brent method in this paper. Generally, there are multiple solutions or roots at a certain frequency or period, and the root corresponding to the lowest c is the fundamental mode while other roots are the higher modes.

3. Computation of PDs of Surface-Wave Phase Velocities

In this section, we derive the algorithms, based on matrix multiplication, for computation of PDs of surface-wave phase velocities. In fact, the surface-wave dispersion equation F also depends on the physical parameters of the subsurface. To facilitate the description, let a vector Q indicate these parameters [35]. In other words,

(2a)$Q_{\mathrm{R}}=\left({\alpha }_1,\cdots ,{\alpha }_{n+1},{\beta }_1,\cdots ,{\beta }_{n+1},{\rho }_1,\cdots ,{\rho }_{n+1},h_1,\cdots ,h_n\right)$

for the case of Rayleigh waves, or

(2b)$Q_{\mathrm{L}}=\left({\beta }_1,\cdots ,{\beta }_{n+1},{\rho }_1,\cdots ,{\rho }_{n+1},h_1,\cdots ,h_n\right)$

for the case of Love waves. Hence, the completed form of F is

(3)$F\left(\omega ,c,Q_{\mathrm{R}}\right)=0\mathrm{or}F\left(\omega ,c,Q_{\mathrm{L}}\right)=0.$

Suppose we want to compute $\partial c/\partial p_t$, where p represents a generic physical parameter such as α, β, ρ, or h of the objective layer t. Making use of the implicit function theory, $\partial c/\partial p_t$ can be written as:

(4)$\partial c/\partial p_t=-\frac{\partial F/\partial p_t}{\partial F/\partial c}.$

In the following, let us first consider $\partial F/\partial c$. On the basis of recursion (1) and the derivative rules, the expression for $\partial F/\partial c$ can be derived as follows:

(5)$\partial F/\partial c=\frac{\partial {\overline{X}}_{n+1}^{*}}{\partial c}{\overline{V}}_{n+1}^{*}+{\overline{X}}_{n+1}^{*}\frac{\partial {\overline{V}}_{n+1}^{*}}{\partial c}.$

Obviously, the calculation of $\partial {\overline{V}}_{n+1}^{*}/\partial c$ is straightforward, whereas $\partial {\overline{X}}_{n+1}^{*}/\partial c$ has to be calculated recursively starting from the first layer via formula (1):

(6)$\partial {\overline{X}}_{m+1}^{*}/\partial c=\frac{\partial {\overline{X}}_m^{*}}{\partial c}{\overline{T}}_m^{*}+{\overline{X}}_m^{*}\frac{\partial {\overline{T}}_m^{*}}{\partial c}, m=1,2,\cdots ,n.$

For a clearer description of the algorithm (Table 1), let us introduce a new vector, i.e., ${\overline{Y}}_m^{*}=\partial {\overline{X}}_m^{*}/\partial c$, thus (6) can be reformulated as follows:

(7)${\overline{Y}}_{m+1}^{*}={\overline{Y}}_m^{*}{\overline{T}}_m^{*}+{\overline{X}}_m^{*}\frac{\partial {\overline{T}}_m^{*}}{\partial c}, m=1,2,\cdots ,n.$

This is obviously a zero vector as ${\overline{X}}_1^{*}$ is a constant vector. A similar algorithm can be derived for $\partial F/\partial p_t$ (Table 2), except that there are two cases that need to be taken into account according to the sequence number t (whether it equals to n + 1). Therefore, by using these two algorithms and (4) in conjunction, the PDs can be computed analytically. We only show the elements of $\partial {\overline{V}}_{n+1}^{*}/\partial c$, $\partial {\overline{T}}_m^{*}/\partial c$, $\partial {\overline{V}}_{n+1}^{*}/\partial p_{n+1}$, and $\partial {\overline{T}}_t^{*}/\partial p_t$ for Love waves in Appendix B. As for Rayleigh waves, due to the lengthy formulas and the limited space, we suggest readers go directly to the source code according to the description in Section 5.

4. Computation of PDs of Surface-Wave Group Velocities

The group velocity is the propagation velocity of the energy of a wave packet and is related to the phase velocity by the following expression [31]:

(8)$U=\partial \omega /\partial k=c/\left[1-\left(\omega /c\right)\partial c/\partial \omega \right]$

where k is the angular wavenumber. Therefore, if certain modes of phase-velocity DCs have been acquired, then the corresponding group-velocity DCs can be calculated by (8).

Two different approaches are implemented in the package: (i) a fully analytical method via implicit function theory (for Love waves); (ii) a hemi-analytical method (for Rayleigh waves). In the following, we will explain both approaches, respectively.

4.1. Fully Analytical Method

The first order derivative $\partial c/\partial \omega$ in (8) also can be calculated by using the implicit function theory:

(9)$\partial c/\partial \omega =-\frac{\partial F/\partial \omega }{\partial F/\partial c}.$

Although the pseudo-code in Table 1 is for $\partial F/\partial c$, it is apparently applicable to $\partial F/\partial \omega$, provided that c is replaced by ω. Similarly, the PDs of group velocity with respect to certain generic physical parameter p of the objective layer t can be derived from (8) [30,31,32,33,36]:

(10)$\frac{\partial U}{\partial p_t}=\frac{U}{c}\left(2-\frac{U}{c}\right)\frac{\partial c}{\partial p_t}+\omega \frac{U^2}{c^2}\frac{{\partial }^{2}c}{\partial \omega \partial p_t}.$

The second-order mixed partial derivative in (10) can be derived by further differentiating (4) with respect to ω [31]:

(11)$\frac{{\partial }^{2}c}{\partial \omega \partial p_t}=-\left[\frac{{\partial }^{2}F}{\partial \omega \partial p_t}+\frac{{\partial }^{2}F}{\partial \omega \partial c}\frac{\partial c}{\partial p_t}+\left(\frac{{\partial }^{2}F}{\partial c\partial p_t}+\frac{{\partial }^{2}F}{\partial c^2}\frac{\partial c}{\partial p_t}\right)\frac{\partial c}{\partial \omega }\right]/\frac{\partial F}{\partial c}.$

In the following, we take $\frac{{\partial }^{2}F}{\partial c^2}$ as an example to illustrate how to calculate the second-order partial derivatives analytically. As for $\frac{{\partial }^{2}F}{\partial c\partial p_t}$, $\frac{{\partial }^{2}F}{\partial \omega \partial p_t}$, and $\frac{{\partial }^{2}F}{\partial \omega \partial c}$, they can be calculated in the same way.

The computational procedure for $\partial F/\partial c$ (Table 1) can be performed in a recursive manner:

(12)$\begin{array}{l}{\overline{X}}_{m+1}^{*}={\overline{X}}_m^{*}{\overline{T}}_m^{*}, m=1,2,\cdots ,n;\\ {\overline{Y}}_{m+1}^{*}={\overline{Y}}_m^{*}{\overline{T}}_m^{*}+{\overline{X}}_m^{*}\frac{\partial {\overline{T}}_m^{*}}{\partial c}, m=1,2,\cdots ,n;\\ \partial F/\partial c={\overline{Y}}_{n+1}^{*}{\overline{V}}_{n+1}^{*}+{\overline{X}}_{n+1}^{*}\frac{\partial {\overline{V}}_{n+1}^{*}}{\partial c}.\end{array}$

Therefore, we differentiate the last expression above with respect to c and thus obtain

(13)${\partial }^{2}F/\partial c^2=\frac{\partial {\overline{Y}}_{n+1}^{*}}{\partial c}{\overline{V}}_{n+1}^{*}+2{\overline{Y}}_{n+1}^{*}\frac{\partial {\overline{V}}_{n+1}^{*}}{\partial c}+{\overline{X}}_{n+1}^{*}\frac{{\partial }^2{\overline{V}}_{n+1}^{*}}{\partial c^2}.$

Then, in order to obtain $\partial {\overline{Y}}_{n+1}^{*}/\partial c$, we shall resort to the second expression in (12):

(14)$\partial {\overline{Y}}_{m+1}^{*}/\partial c=\frac{\partial {\overline{Y}}_m^{*}}{\partial c}{\overline{T}}_m^{*}+2{\overline{Y}}_m^{*}\frac{\partial {\overline{T}}_m^{*}}{\partial c}+{\overline{X}}_m^{*}\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial c^2},m=1,2,\cdots ,n.$

Similarly, if we introduce a new vector ${\overline{Z}}_m^{*}=\partial {\overline{Y}}_m^{*}/\partial c$, then (14) can be reformulated as follows:

(15)${\overline{Z}}_{m+1}^{*}={\overline{Z}}_m^{*}{\overline{T}}_m^{*}+2{\overline{Y}}_m^{*}\frac{\partial {\overline{T}}_m^{*}}{\partial c}+{\overline{X}}_m^{*}\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial c^2},m=1,2,\cdots ,n.$

Obviously, ${\overline{Z}}_1^{*}$ is a zero vector as well. Thus, by taking advantage of (12), (15), and (13), we can finally implement the analytical computation of ${\partial }^{2}F/\partial c^2$ (Table 3). See Appendix C for the elements of ${\partial }^2{\overline{T}}_m^{*}/\partial c^2$ and of ${\partial }^2{\overline{V}}_{n+1}^{*}/\partial c^2$ for Love waves.

4.2. Hemi-Analytical Method

In the fully analytical method, the crucial second-order partial derivative term $\frac{{\partial }^{2}c}{\partial \omega \partial p_t}$ in (10) is substituted by (11) and then solved via some recursive procedures. In contrast, the hemi-analytical method approximates this term with a central finite-difference scheme [36], a procedure which will be briefly reviewed in the following.

Suppose that the partial derivative $\partial U/\partial p_t$ is desired at the angular frequency ${\omega }_0$ and that the phase velocity $c_0$ and group velocity $U_0$ are computed at ${\omega }_0$. Define two other frequencies around ${\omega }_0$: ${\omega }_{+1}={\omega }_0{e}^{+\delta }$ and ${\omega }_{-1}={\omega }_0{e}^{-\delta }$. $\delta$ is a small quantity, denoting a perturbation of log $\omega$; see [36] for a more detailed explanation. And suppose the quantities $c_{+1}$, $U_{+1}$, and $\partial c_{+1}/\partial p_t$ are corresponding to ${\omega }_{+1}$ and have been obtained. The same hypothesis is applied to the quantities $c_{-1}$, $U_{-1}$, and $\partial c_{-1}/\partial p_t$ at ${\omega }_{-1}$. Then the PD of the group velocity at ${\omega }_0$ can be computed via

(16)$\partial U_0/\partial p_t=\frac{U_0}{c_0}\left(2-\frac{U_0}{c_0}\right)\frac{\partial c_0}{\partial p_t}+\frac{U_0^2}{c_0^2}\left(\frac{\partial c_{+1}/\partial p_t-\partial c_{-1}/\partial p_t}{2\delta }\right)$

We call (16) ‘hemi-analytical’ due to the fact that all the first-order partial derivatives and $U_0$ in (16) are computed by the analytical approach discussed previously, leaving only $\frac{{\partial }^{2}c}{\partial \omega \partial p_t}$ approximated by the central finite-difference.

5. Package Description

SWPD, written in MATLAB, consists of two independent parts for dealing with Rayleigh and Love wave respectively. Most source files of the package are written for implementing the partial differentiations of matrix ${\overline{T}}_m^{*}$ and of dispersion function F. To make it clearer, we list these files in Table A1 and Table A2, respectively. As shown in these two Tables, for example, reduced_delta_love.m and reduced_delta.m implement the delta matrix recursion (1) for Love and Rayleigh waves, respectively. Files reduced_delta_love_dfdc.m and reduced_delta_dfdc.m perform the pseudo-code in Table 1 to yield the PDs of the dispersion function F with respect to c. When $\partial F/\partial \beta$ is desired (as the pseudo-code shown in Table 2), one can go to files reduced_delta_love_dfdb.m and reduced_delta_dfdb.m. Similarly, reduced_delta_love_dfdcc.m has implemented the algorithm in Table 3. In summary, all the main *.m files that perform partial differentiations are summarized in Table A2 and they need to call the sub-function files in Table A1 to accomplish the relevant computations.

Besides those scripts described above, there are a few other *.m files deserving to be mentioned. rayleighphase.m and lovephase.m are used to model the multimodal phase-velocity DCs. The hemi-analytical method has been implemented in four *.m files whose filenames contain the string ‘rodi’. In fact, we have included five demo scripts, whose filenames begin with the string ‘demo’. Hence, one can replicate all figures of this paper by running the demo scripts directly and so to become familiar with SWPD as quickly as possible.

6. Examples

6.1. Love-Wave Case

First, a two-layer model (Model Nov71 in Table 4) is chosen to demonstrate the validity of the proposed algorithms and of the developed code for the case of Love waves. Figure 2 shows the computational results for Model Nov71. Although the package is capable of modeling multiple modes, we only show the results of the fundamental mode for comparison with the published studies [30,31]. It is noteworthy that, in Figure 2 we plot curves of $10\times \partial c/\partial h$ and $10\times \partial U/\partial h$ instead of their original values due to their large differences in quantitative values compared to the other derivatives, so to make the changes of these two curves recognizable. In addition, in the work of Novotný [31], to facilitate the derivation of the analytical formulism specifically for the case of two layers, he introduced a dimensionless quantity $\rho ={\rho }_2/{\rho }_1$ and plotted the curves of $\partial c/\partial \rho$ and $\partial U/\partial \rho$. Therefore, to obtain the same curves, we have made use of the chain rule:

(17)$\partial c/\partial \rho =\frac{\partial c}{\partial {\rho }_2}\frac{\partial {\rho }_2}{\partial \rho }={\rho }_1\frac{\partial c}{\partial {\rho }_2} \mathrm{and} \partial U/\partial \rho =\frac{\partial U}{\partial {\rho }_2}\frac{\partial {\rho }_2}{\partial \rho }={\rho }_1\frac{\partial U}{\partial {\rho }_2}.$

Obviously, the computational results in Figure 2 are rather consistent with the work of Novotný [31] (see Figure 1, Figure 2 and Figure 3 of his paper). To verify the correctness of these results more clearly, the relevant quantities calculated at 20, 30, and 40 s are shown in Table 5 (related to phase velocity) and in Table 6 (related to group velocity), respectively. As shown in the Tables, they are in good agreement with the outcomes from Novotný [31] (see Table 1 and Table 2 of his paper).

Next, an eight-layer model (Model Nov70 in Table 4) is used for a further verification. The PDs of the fundamental-mode phase and group velocity are presented in Figure 3. Also, the results calculated at 20 and 40 s for individual layer are listed in Table 7 and Table 8 for detailed comparisons. It should be pointed out that Novotný [30] (Figure 1, Figure 2, Figure 3, Figure 4 and Figure 5 in that paper) did not plot some partial derivatives because of their too small absolute values, whereas we plot all partial derivatives, regardless of their magnitude. After comparison, we can conclude that the results for the case of multiple layers are also satisfactory.

6.2. Rayleigh-Wave Case

6.2.1. A Two-Layer Model

Let us further consider Rayleigh waves. In the first example, we adopt a two-layer model that is almost the same with Model Nov71, except that the P-wave velocities of the layers are also taken into account (${\alpha }_1$ = 6000 m∙s⁻¹ and ${\alpha }_2$ = 8000 m∙s⁻¹) [32]. Figure 4 displays the related computational outcomes of fundamental-mode phase and group velocity. They show apparently good agreements with [32] (Figure 1, Figure 2 and Figure 3 in that paper). Table 9 also shows the relevant results of phase velocity at three different periods for a detailed comparison. It is worth to point out that although we have used difference approximation while coping with the PDs of group velocity, the quantities acquired by the hemi-analytical approach still can be accurate to five decimal places as compared to the fully analytical approach [32] (see Figure 4c and Table 10). Hence, the hemi-analytical approach is a rather effective and accurate scheme that can be used to construct the Jacobian matrix while the linearized inversion of Rayleigh-wave group-velocity DCs is to be conducted.

6.2.2. A Six-Layer Model

All of the above-mentioned models are aimed at the scale of crustal seismology. In this example, we chose a near-surface six-layer model (Table 11) to demonstrate the applicability and effectiveness of our code at the scale of near-surface geophysics. This model was also chosen because some researchers have reported their results [1,35].

As shown in Figure 5, the analytical PDs of the fundamental-mode Rayleigh-wave phase velocity of the six-layer model are quite consistent with [35] (Figure 3, Figure 4 and Figure 5 of his paper). In order to show more details, as an example, the PDs with respect to the S-wave velocities at five different frequencies, along with the outcomes from other researchers, are presented in Table 12 (modified from [35]). Individual column of these matrices corresponds to a layer, whereas individual row corresponds to a frequency. Obviously, only the numerical results from Lai and Rix [6] show a certain lack of accuracy due to the adopted numerical integration in their code [35], whereas other results (including our result) are quite similar.

7. Discussion and Conclusions

The accurate and reliable calculation of PDs of surface-wave DCs is an essential component of any linearized inversion strategy based on surface-wave DCs. In this paper, we have rederived a series of algorithms and the related expressions based on the reduced delta matrix theory to compute analytically the PDs of surface-wave (both Love and Rayleigh waves) phase-velocity DCs. Moreover, we have also applied this analytical method to compute the PDs of Love-wave group velocity. Although Equations (9)–(15) constitute an abstract framework applicable to both waves, the derivation of all the second-order mixed PDs for Rayleigh waves could be rather difficult. Therefore, we have implemented a high-precision hemi-analytical method, which computes all the first-order partial differentiation analytically (Section 3) and approximates the only mixed second-order partial differentiation term with a central finite-difference scheme. For the fully analytical computation of Rayleigh-wave group-velocity PDs, Novotný et al. [32] is the only work we could find, despite the fact that only the case of two-layer model was studied due to the complexity of this problem. The comparison (Section 6.2.1) demonstrates that this mixed computational strategy is satisfying, with a comparable accuracy to five decimal places, as shown in Table 10.

An open-source MATLAB package named SWPD has been developed accompanying this paper. All examples in this paper have been validated by comparing with the results of previous studies, and all the results computed by our package are satisfactory with rather high precision. Therefore, it has been demonstrated that the analytical expressions and the algorithms derived in this paper are effective. Also, these examples can be reproduced by other users with the help of four demo scripts in the package. Upon becoming familiar with this package, people could modify and extend it without any difficulty, including any further development of the surface-wave linearized inversion strategy. We think this package could be a useful tool for people who are engaged in the field of near-surface geophysics, especially for those who are interested in surface-wave exploration methods.

Author Contributions

Conceptualization, D.W., X.W. and Q.S.; Methodology, D.W., X.W. and Q.S.; Software, D.W.; Formal analysis, D.W., X.W., Q.S. and T.Z.; Funding acquisition, X.W. and Q.S.; Investigation, D.W., Q.S. and T.Z.; Project administration, X.W., Q.S. and T.Z.; Resources, X.W. and Q.S.; Validation, D.W., X.W., Q.S. and T.Z.; Writing—original draft, D.W.; Writing—review & editing D.W., X.W., Q.S. and T.Z.

Funding

This research is funded by the National Science and Technology Major Project of the Ministry of Science and Technology of China “the study of key geophysical exploration technologies for Lower Palaeozoic Erathem-Precambrian basement tectonics (2016ZX05004-003)”

Acknowledgments

We are grateful to the valuable comments by the anonymous reviewers, which have helped improve this manuscript.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

The expressions presented below are adapted from Buchen & Ben-Hador [23] just for convenience of reference.

Basic physical quantities employed:

• $\omega =$ angular frequency

• $c=$ phase velocity

• $k=\omega /c=$ angular wavenumber

• $i=\sqrt{-1}$ = imaginary unit

Basic derived quantities employed for the m-th layer (subscript m is omitted for simplicity):

(A1)

$r=\{\begin{array}{ll}\sqrt{1-c^2/{\alpha }^2},& c<\alpha \\ i\sqrt{c^2/{\alpha }^2-1},& c>\alpha \end{array}$

(A2)

$r=\{\begin{array}{ll}\sqrt{1-c^2/{\alpha }^2},& c<\alpha \\ i\sqrt{c^2/{\alpha }^2-1},& c>\alpha \end{array}$

(A3)

$C_{\alpha }=\{\begin{array}{ll}\cosh \left(krh\right),& c<\alpha \\ \cos \left(ikrh\right),& c>\alpha \end{array}$

(A4)

$C_{\beta }=\{\begin{array}{ll}\cosh \left(ksh\right),& c<\beta \\ \cos \left(iksh\right),& c>\beta \end{array}$

(A5)

$S_{\alpha }=\{\begin{array}{ll}\sinh \left(krh\right),& c<\alpha \\ -i\sin \left(ikrh\right),& c>\alpha \end{array}$

(A6)

$S_{\beta }=\{\begin{array}{ll}\sinh \left(ksh\right),& c<\beta \\ -i\sin \left(iksh\right),& c>\beta \end{array}$

(A7)

$\gamma ={\beta }^2/c^2, t=2-1/\gamma , \mu =\rho {\beta }^2$

(A8)

$\begin{array}{ll}ForLovewaves:& ForRayleighwaves:\\ {\overline{X}}_1^{*}=\left[\begin{array}{cc}0& 1\end{array}\right]& {\overline{X}}_1^{*}=\left[\begin{array}{ccccc}0& 0& 0& 0& 1\end{array}\right]\\ {\overline{T}}_m^{*}={\left[\begin{array}{cc}{C}_{\beta }& \frac{S_{\beta }}{\mu s}\\ \mu s{S}_{\beta }& C_{\beta }\end{array}\right]}_m& {\overline{T}}_m^{*}={\left[\begin{array}{ccccc}{\overline{T}}_{11}& {\overline{T}}_{12}& {\overline{T}}_{13}& {\overline{T}}_{14}& {\overline{T}}_{16}\\ 2{\overline{T}}_{21}& 2{\overline{T}}_{22}-1& 2{\overline{T}}_{23}& 2{\overline{T}}_{24}& 2{\overline{T}}_{26}\\ {\overline{T}}_{31}& {\overline{T}}_{32}& {\overline{T}}_{33}& {\overline{T}}_{34}& {\overline{T}}_{36}\\ {\overline{T}}_{41}& {\overline{T}}_{42}& {\overline{T}}_{43}& {\overline{T}}_{44}& {\overline{T}}_{46}\\ {\overline{T}}_{61}& {\overline{T}}_{62}& {\overline{T}}_{63}& {\overline{T}}_{64}& {\overline{T}}_{66}\end{array}\right]}_m\\ {\overline{V}}_{n+1}^{*}={\left[\begin{array}{cc}1& \mu s\end{array}\right]}_{n+1}^T& {\overline{V}}_{n+1}^{*}={\left[\begin{array}{ccccc}1-rs& 2\mu \left(t-2rs\right)& \mu s\left(2-t\right)& -\mu r\left(2-t\right)& {\mu }^2\left(4rs-t^2\right)\end{array}\right]}_{n+1}^T\end{array}$

where the superscript T means the matrix transpose.

The elements of matrix ${\overline{T}}_m^{*}$ for Rayleigh wave:

(A9)

$Q_j=\left(\frac{t^j}{rs}+2^{j}rs\right)S_{\alpha }{S}_{\beta }, 0\le j\le 4$

(A10)

${\overline{T}}_{11}={\gamma }^2\left[-4t+\left(t^2+4\right)C_{\alpha }{C}_{\beta }-Q_2\right]$

(A11)

${\overline{T}}_{12}=\frac{{\gamma }^2}{\mu }\left[\left(2+t\right)\left(1-C_{\alpha }{C}_{\beta }\right)+Q_1\right]$

(A12)

${\overline{T}}_{13}=\frac{\gamma }{\mu }\left(C_{\alpha }{S}_{\beta }/s-r{S}_{\alpha }{C}_{\beta }\right)$

(A13)

${\overline{T}}_{14}=\frac{\gamma }{\mu }\left(s{C}_{\alpha }{S}_{\beta }-S_{\alpha }{C}_{\beta }/r\right)$

(A14)

${\overline{T}}_{16}=\frac{{\gamma }^2}{{\mu }^2}\left[2\left(1-C_{\alpha }{C}_{\beta }\right)+Q_0\right]$

(A15)

${\overline{T}}_{21}=\mu {\gamma }^2\left[-2t\left(t+2\right)\left(1-C_{\alpha }{C}_{\beta }\right)-Q_3\right]$

(A16)

${\overline{T}}_{22}=1+C_{\alpha }{C}_{\beta }-{\overline{T}}_{11}$

(A17)

${\overline{T}}_{23}=\gamma \left(t{C}_{\alpha }{S}_{\beta }/s-2r{S}_{\alpha }{C}_{\beta }\right)$

(A18)

${\overline{T}}_{24}=\gamma \left(2s{C}_{\alpha }{S}_{\beta }-t{S}_{\alpha }{C}_{\beta }/r\right)$

(A19)

${\overline{T}}_{26}={\overline{T}}_{12}$

(A20)

${\overline{T}}_{31}=\mu \gamma \left(4s{C}_{\alpha }{S}_{\beta }-t^2{S}_{\alpha }{C}_{\beta }/r\right)$

(A21)

${\overline{T}}_{32}=-{\overline{T}}_{24}$

(A22)

${\overline{T}}_{33}=C_{\alpha }{C}_{\beta }$

(A23)

${\overline{T}}_{34}=-s{S}_{\alpha }{S}_{\beta }/r$

(A24)

${\overline{T}}_{36}=-{\overline{T}}_{14}$

(A25)

${\overline{T}}_{41}=\mu \gamma \left(t^2{C}_{\alpha }{S}_{\beta }/s-4r{S}_{\alpha }{C}_{\beta }\right)$

(A26)

${\overline{T}}_{42}=-{\overline{T}}_{23}$

(A27)

${\overline{T}}_{43}=-r{S}_{\alpha }{S}_{\beta }/s$

(A28)

${\overline{T}}_{44}={\overline{T}}_{33}$

(A29)

${\overline{T}}_{46}=-{\overline{T}}_{13}$

(A30)

${\overline{T}}_{61}={\mu }^2{\gamma }^2\left[8t^2\left(1-C_{\alpha }{C}_{\beta }\right)+Q_4\right]$

(A31)

${\overline{T}}_{62}={\overline{T}}_{21}$

(A32)

${\overline{T}}_{63}=-{\overline{T}}_{41}$

(A33)

${\overline{T}}_{64}=-{\overline{T}}_{31}$

(A34)

${\overline{T}}_{66}={\overline{T}}_{11}$

Appendix B

The expressions listed below are used for the construction of the first-order PDs of matrices ${\overline{T}}_m^{*}$ and ${\overline{V}}_{n+1}^{*}$ for Love waves. For clarity, we use the same notation as Cercato [35], i.e., the prime mark (e.g., $s^{\prime }$) is used exclusively to indicate PD with respect to phase velocity c, whereas any subscript indicates PD with respect to that parameter (e.g., $s_{\beta }=\partial s/\partial \beta$ and ${\left(C_{\beta }\right)}_{\beta }=\partial C_{\beta }/\partial \beta$).

Define an intermediate variable

(A35)

$s_0=1-c^2/{\beta }^2$

For $\partial {\overline{T}}_m^{*}/\partial c$:

(A36)

$s_0^{\prime }=-2c/{\beta }^2$

(A37)

$s^{\prime }=0.5s_0^{\prime }/s$

(A38)

$C_{\beta }^{\prime }=S_{\beta }kh\left(s^{\prime }-s/c\right)$

(A39)

$s_{\beta }^{\prime }=C_{\beta }kh\left(s^{\prime }-s/c\right)$

(A40)

$\partial {\overline{T}}_m^{*}/\partial c=\left[\begin{array}{cc}{C}_{\beta }^{\prime }& \frac{S_{\beta }^{\prime }s-S_{\beta }{s}^{\prime }}{\mu s^2}\\ \mu \left(s^{\prime }{S}_{\beta }+s{s}_{\beta }^{\prime }\right)& C_{\beta }^{\prime }\end{array}\right]$

For $\partial {\overline{V}}_{n+1}^{*}/\partial c$:

(A41)

$\partial {\overline{V}}_{n+1}^{*}/\partial c={\left[\begin{array}{cc}0& s_{\mu }^{\prime }\end{array}\right]}^T$

For $\partial {\overline{T}}_m^{*}/\partial \beta$:

(A42)

${\left(s_0\right)}_{\beta }=2c^2/{\beta }^3$

(A43)

${\mu }_{\beta }=2\rho \beta$

(A44)

$s_{\beta }=0.5{\left(s_0\right)}_{\beta }/s$

(A45)

${\left(C_{\beta }\right)}_{\beta }=kh{S}_{\beta }{s}_{\beta }, {\left(S_{\beta }\right)}_{\beta }=kh{C}_{\beta }{s}_{\beta }$

(A46)

$\partial {\overline{T}}_m^{*}/\partial \beta =\left[\begin{array}{cc}{\left(C_{\beta }\right)}_{\beta }& -\frac{{\mu }_{\beta }{S}_{\beta }}{{\mu }^{2}s}+\frac{{\left(S_{\beta }\right)}_{\beta }s-S_{\beta }{s}_{\beta }}{\mu s^2}\\ {\mu }_{\beta }s{S}_{\beta }+\mu \left(S_{\beta }{s}_{\beta }+s{\left(S_{\beta }\right)}_{\beta }\right)& {\left(C_{\beta }\right)}_{\beta }\end{array}\right]$

For $\partial {\overline{V}}_{n+1}^{*}/\partial \beta$:

(A48)

$\partial {\overline{V}}_{n+1}^{*}/\partial \beta ={\left[\begin{array}{cc}0& {\mu }_{\beta }s+\mu s_{\beta }\end{array}\right]}^T$

For $\partial {\overline{T}}_m^{*}/\partial \rho$:

(A49)

${\mu }_{\rho }={\beta }^2$

(A50)

$\partial {\overline{T}}_m^{*}/\partial \rho =\left[\begin{array}{cc}0& -\frac{{\mu }_{\rho }{S}_{\beta }}{{\mu }^{2}s}\\ {\mu }_{\rho }s{S}_{\beta }& 0\end{array}\right]$

For $\partial {\overline{V}}_{n+1}^{*}/\partial \rho$:

(A51)

$\partial {\overline{V}}_{n+1}^{*}/\partial \rho ={\left[\begin{array}{cc}0& {\mu }_{\rho }s\end{array}\right]}^T$

For $\partial {\overline{T}}_m^{*}/\partial h$:

(A52)

${\left(C_{\beta }\right)}_h=S_{\beta }ks, {\left(S_{\beta }\right)}_h=C_{\beta }ks$

(A53)

$\partial {\overline{T}}_m^{*}/\partial h=\left[\begin{array}{cc}{\left(C_{\beta }\right)}_h& \frac{{\left(S_{\beta }\right)}_h}{\mu s}\\ \mu s{\left(S_{\beta }\right)}_h& {\left(C_{\beta }\right)}_h\end{array}\right]$

Appendix C

The double-prime mark is used for indicating the second-order PDs with respect to c for Love waves.

(A54)

$s_0^{″}=-2/{\beta }^2$

(A55)

$s^{″}=\left(0.5s_0^{″}-{\left(s^{\prime }\right)}^2\right)/s$

(A56)

${C^{″}}_{\beta }=kh\left\{\left({S^{\prime }}_{\beta }-\frac{S_{\beta }}{c}\right)\left(s^{\prime }-\frac{s}{c}\right)+S_{\beta }\left[s^{″}-\left(s^{\prime }-\frac{s}{c}\right)\frac{1}{c}\right]\right\}$

(A57)

${S^{″}}_{\beta }=kh\left\{\left({C^{\prime }}_{\beta }-\frac{C_{\beta }}{c}\right)\left(s^{\prime }-\frac{s}{c}\right)+C_{\beta }\left[s^{″}-\left(s^{\prime }-\frac{s}{c}\right)\frac{1}{c}\right]\right\}$

For ${\partial }^2{\overline{T}}_m^{*}/\partial c^2$:

(A58)

$y={S^{\prime }}_{\beta }s-S_{\beta }{s}^{\prime }$

(A59)

$y^{\prime }={S^{″}}_{\beta }s-S_{\beta }{s}^{″}$

(A60)

${\partial }^2{\overline{T}}_m^{*}/\partial c^2=\left[\begin{array}{cc}{C^{″}}_{\beta }& \frac{y^{\prime }{s}^2-2ys{s}^{\prime }}{\mu s^4}\\ \mu \left(s^{″}{S}_{\beta }+2s^{\prime }{S^{\prime }}_{\beta }+s{S^{″}}_{\beta }\right)& {C^{″}}_{\beta }\end{array}\right]$

For ${\partial }^2{\overline{V}}_{n+1}^{*}/\partial c^2$:

(A61)

${\partial }^2{\overline{V}}_{n+1}^{*}/\partial c^2={\left[\begin{array}{cc}0& s^{″}\mu \end{array}\right]}^T$

Appendix D

Table A1

MATLAB *.m files that compute matrix ${\overline{T}}_m^{*}$ and its partial derivatives (PDs).

Role	Love waves	Rayleigh waves
${\overline{T}}_m^{*}$	get_T_love.m	get_T.m
$\partial {\overline{T}}_m^{*}/\partial \alpha$		get_dTda.m
$\partial {\overline{T}}_m^{*}/\partial \beta$	get_dTdb_love.m	get_dTdb.m
$\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial \beta \partial c}$	get_dTdbc_love.m
$\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial \beta \partial \omega }$	get_dTdbw_love.m
$\partial {\overline{T}}_m^{*}/\partial c$	get_dTdc_love.m	get_dTdc.m
${\partial }^2{\overline{T}}_m^{*}/\partial c^2$	get_dTdcc_love.m
$\partial {\overline{T}}_m^{*}/\partial \omega$	get_dTdw_love.m	get_dTdw.m
$\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial \omega \partial c}$	get_dTdwc_love.m
$\partial {\overline{T}}_m^{*}/\partial h$	get_dTdh_love.m	get_dTdh.m
$\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial h\partial c}$	get_dTdhc_love.m
$\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial h\partial \omega }$	get_dTdhw_love.m
$\partial {\overline{T}}_m^{*}/\partial \rho$	get_dTdrho_love.m	get_dTdrho.m
$\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial \rho \partial c}$	get_dTdrhoc_love.m
$\frac{{\partial }^2{\overline{T}}_m^{*}}{\partial \rho \partial \omega }$	get_dTdrhow_love.m

Table A2

MATLAB *.m files that compute dispersion function F and its PDs.

Role	Love waves	Rayleigh waves
$F$	reduced_delta_love.m	reduced_delta.m
$\partial F/\partial \alpha$		reduced_delta_dfda.m
$\partial F/\partial \beta$	reduced_delta_love_dfdb.m	reduced_delta_dfdb.m
$\frac{{\partial }^{2}F}{\partial \beta \partial c}$	reduced_delta_love_dfdbc.m
$\frac{{\partial }^{2}F}{\partial \beta \partial \omega }$	reduced_delta_love_dfdbw.m
$\partial F/\partial c$	reduced_delta_love_dfdc.m	reduced_delta _dfdc.m
${\partial }^{2}F/\partial c^2$	reduced_delta_love_dfdcc.m
$\partial F/\partial h$	reduced_delta_love_dfdh.m	reduced_delta _dfdh.m
$\frac{{\partial }^{2}F}{\partial h\partial c}$	reduced_delta_love_dfdhc.m
$\frac{{\partial }^{2}F}{\partial h\partial \omega }$	reduced_delta_love_dfdhw.m
$\partial F/\partial \rho$	reduced_delta_love_dfdrho.m	reduced_delta _dfdrho.m
$\frac{{\partial }^{2}F}{\partial \rho \partial c}$	reduced_delta_love_dfdrhoc.m
$\frac{{\partial }^{2}F}{\partial \rho \partial \omega }$	reduced_delta_love_dfdrhow.m
$\partial F/\partial \omega$	reduced_delta_love_dfdw.m	reduced_delta _dfdw.m
$\frac{{\partial }^{2}F}{\partial \omega \partial c}$	reduced_delta_love_dfdwc.m

Figures and Tables

Figure 1. A multilayered earth model with each layer being homogeneous, isotropic and elastic.

Figure 2. Computational results of Model Nov71: (a) phase and group velocities of the fundamental mode; (b) Partial derivatives (PDs) of the phase velocity; (c) PDs of the group velocity.

Figure 2. Computational results of Model Nov71: (a) phase and group velocities of the fundamental mode; (b) Partial derivatives (PDs) of the phase velocity; (c) PDs of the group velocity.

Figure 3. Computational results of Model Nov70. PDs of the phase velocity with respect to β (a), ρ (c), and h (e) are presented in the left column, whereas (b,d,f) show the corresponding PDs of the group velocity. Each line represents the PD of that numbered layer.

Figure 4. Computational results of Model Nov71 for Rayleigh wave: (a) phase and group velocities for the fundamental mode; (b) analytical PDs of the fundamental-mode phase velocity; (c) the PDs of the fundamental-mode group velocity acquired by the hemi-analytical approach.

Figure 5. Analytical PDs of the fundamental-mode phase velocity with respect to β (a), α (b), ρ (c), and h (d) for the six-layer model in Table 10.

Table 1

Pseudo-code for computation of $\partial F/\partial c$.

Pseudo-code 1 The computation of $\partial F/\partial c$
Input: angular frequency ω, phase velocity c, earth model $Q_{\mathrm{R}}$ or $Q_{\mathrm{L}}$
1:	Initialize four vectors: xold = ${\overline{X}}_1^{*}$, yold = ${\overline{Y}}_1^{*}$, xnew = 0, and ynew = 0
2:	form = 1–n do
3:	calculate ${\overline{T}}_m^{*}$ and $\partial {\overline{T}}_m^{*}/\partial c$
4:	update xnew: xnew = xold $\times {\overline{T}}_m^{*}$
5:	update ynew: ynew = yold $\times {\overline{T}}_m^{*}$ + xold $\times \partial {\overline{T}}_m^{*}/\partial c$
6:	update xold: xold = xnew
7:	update yold: yold = ynew
8:	end for
9:	compute ${\overline{V}}_{n+1}^{*}$ and $\partial {\overline{V}}_{n+1}^{*}/\partial c$
Output:$\partial F/\partial c$ = ynew $\times {\overline{V}}_{n+1}^{*}$ + xnew $\times \partial {\overline{V}}_{n+1}^{*}/\partial c$

Table 2

Pseudo-code for computation of $\partial F/\partial p_t$.

Pseudo-code 2 The computation of $\partial F/\partial p_t$
Input: angular frequency ω, phase velocity c, earth model $Q_{\mathrm{R}}$ or $Q_{\mathrm{L}}$, and the sequence number t of the objective layer
1:	Initialize one vector and one matrix: x = ${\overline{X}}_1^{*}$, T = 0
2:	ift equals to n + 1 then
3:	for m = 1–n do
4:	calculate ${\overline{T}}_m^{*}$: T = ${\overline{T}}_m^{*}$
5:	update x: x = x × T
6:	end for
7:	Calculate $\partial {\overline{V}}_{n+1}^{*}/\partial p_{n+1}$
Output:$\partial F/\partial p_{n+1}$ = x $\times \partial {\overline{V}}_{n+1}^{*}/\partial p_{n+1}$
8:	else
9:	for m = 1–n do
10:	if m equals to t then
11:	Calculate $\partial {\overline{T}}_t^{*}/\partial p_t$: T = $\partial {\overline{T}}_t^{*}/\partial p_t$
12:	else
13:	Calculate ${\overline{T}}_m^{*}$: T = ${\overline{T}}_m^{*}$
14:	end if
15:	update x: x = x × T
16:	end for
17:	compute ${\overline{V}}_{n+1}^{*}$
Output:$\partial F/\partial p_t$ = x $\times {\overline{V}}_{n+1}^{*}$
18:	end if

Table 3

Pseudo-code for computation of ${\partial }^{2}F/\partial c^2$.

Pseudo-code 3 The computation of ${\partial }^{2}F/\partial c^2$
Input: angular frequency ω, phase velocity c, earth model $Q_{\mathrm{R}}$ or $Q_{\mathrm{L}}$
1:	Initialize six vectors: xold = ${\overline{X}}_1^{*}$, yold = ${\overline{Y}}_1^{*}$, zold = ${\overline{Z}}_m^{*}$, xnew = 0, ynew = 0, and znew = 0
2:	form = 1–n do
3:	calculate ${\overline{T}}_m^{*}$, $\partial {\overline{T}}_m^{*}/\partial c$ and ${\partial }^2{\overline{T}}_m^{*}/\partial c^2$
4:	update xnew: xnew = xold $\times {\overline{T}}_m^{*}$
5:	update ynew: ynew = yold $\times {\overline{T}}_m^{*}$ + xold $\times \partial {\overline{T}}_m^{*}/\partial c$
6:	update znew: znew = zold $\times {\overline{T}}_m^{*}+2\times \mathrm{yold}\times \partial {\overline{T}}_m^{*}/\partial c+xold\times {\partial }^2{\overline{T}}_m^{*}/\partial c^2$
7:	update xold: xold = xnew
8:	update yold: yold = ynew
9:	update zold: zold = znew
10:	end for
11:	compute ${\overline{V}}_{n+1}^{*}$, $\partial {\overline{V}}_{n+1}^{*}/\partial c$ and ${\partial }^2{\overline{V}}_{n+1}^{*}/\partial c^2$
Output:${\partial }^{2}F/\partial c^2$ = znew $\times {\overline{V}}_{n+1}^{*}+2\times ynew\times \partial {\overline{V}}_{n+1}^{*}/\partial c+xnew\times {\partial }^2{\overline{V}}_{n+1}^{*}/\partial c^2$

Table 4

Two models chosen for validation of the package.

Layer	S-Wave Velocity (m∙s−1)	Density (kg∙m−3)	Thickness (m)
Model Nov71
Two layered model in Novotný (1971) [31]
1	3500	2700	35,000
2	4500	3300	∞
Model Nov70
The model of the Canadian shield CANSD in Novotný (1970) [30]
1	3470	2700	6000
2	3640	2800	10,500
3	3850	2850	18,700
4	4720	3300	80,000
5	4540	3440	100,000
6	4510	3530	100,000
7	4760	3600	80,000
8	5120	3760	∞

Table 5

Computational results relevant to phase velocities for Model Nov71.

T (s)	c (m∙s−1)		∂c/∂β1		∂c/∂β2
	Novotný	This paper	Novotný	This paper	Novotný	This paper
20	3790.45	3790.4529	1.05344	1.0534372	0.12405	0.1240473
30	4010.57	4010.5735	0.88083	0.8808288	0.33920	0.3392031
40	4176.65	4176.6474	0.62679	0.6267904	0.55832	0.5583233
			∂c/∂h		∂c/∂ρ
			Novotný	This paper	Novotný	This paper
20			−0.01299	−0.0129940	0.10281	0.1028084
30			−0.01711	−0.0171069	0.21306	0.2130622
40			−0.01513	−0.0151307	0.25015	0.2501495

Table 6

Computational results relevant to group velocities for Model Nov71.

T (s)	U (m∙s−1)		∂U/∂β1		∂U/∂β2
	Novotný	This paper	Novotný	This paper	Novotný	This paper
20	3384.38	3384.3839	1.16642	1.1664245	−0.14784	−0.1478377
30	3489.61	3489.6075	1.42265	1.4226530	−0.20843	−0.2084253
40	3706.66	3706.6645	1.36767	1.3676722	−0.03879	−0.0387899
			∂U/∂h		∂U/∂ρ
			Novotný	This paper	Novotný	This paper
20			−0.00094	−0.0009381	−0.08530	−0.0853015
30			−0.01576	−0.0157647	0.02334	0.0233376
40			−0.02588	−0.0258753	0.25547	0.2554651

Table 7

Computational results relevant to phase velocities for Model Nov70.

T (s)	∂c/∂βm		∂c/∂ρm		∂c/∂hm
	Novotný	This Paper	Novotný	This Paper	Novotný	This Paper
20	0.25048	0.2504840	−0.05192	−0.0519198	−0.02455	−0.0245543
	0.40458	0.4045771	−0.03721	−0.0372102	−0.01807	−0.0180673
	0.40666	0.4066578	0.03642	0.0364167	−0.01113	−0.0111329
	0.13636	0.1363630	0.04255	0.0425454	0.00000	0.0000024
	0.00022	0.0002246	0.00005	0.0000537	0.00000	0.0000000
	0.00000	0.0000001	0.00000	0.0000000	0.00000	0.0000000
	0.00000	0.0000000	0.00000	0.0000000	0.00000	0.0000000
	0.00000	0.0000000	0.00000	0.0000000
40	0.09234	0.0923393	−0.03562	−0.0356177	−0.01633	−0.0163304
	0.16802	0.1680152	−0.04359	−0.0435873	−0.01415	−0.0141500
	0.26171	0.2617060	−0.02070	−0.0206981	−0.01100	−0.0109961
	0.48998	0.4899821	0.07374	0.0737381	0.00037	0.0003742
	0.10746	0.1074640	0.00813	0.0081273	0.00000	−0.0000047
	0.01742	0.0174226	0.00136	0.0013618	−0.00001	−0.0000138
	0.00147	0.0014714	0.00027	0.0002728	0.00000	−0.0000016
	0.00011	0.0001118	0.00003	0.0000315

Table 8

Computational results relevant to group velocities for Model Nov70.

T (s)	∂U/∂βm		∂U/∂ρm		∂U/∂hm
	Novotný	This Paper	Novotný	This Paper	Novotný	This Paper
20	0.39610	0.3960951	−0.03708	−0.0370822	−0.01867	−0.0186684
	0.56217	0.5621708	0.01013	0.0101325	−0.00798	−0.0079819
	0.29970	0.2996970	0.06946	0.0694558	0.00043	0.0004303
	−0.17980	−0.1797986	−0.03776	−0.0377645	−0.00002	−0.0000189
	−0.00205	−0.0020513	−0.00046	−0.0004573	0.00000	0.0000000
	0.00000	−0.0000027	0.00000	−0.0000006	0.00000	0.0000000
	0.00000	0.0000000	0.00000	0.0000000	0.00000	0.0000000
	0.00000	0.0000000	0.00000	0.0000000
40	0.26627	0.2662661	−0.08695	−0.0869520	−0.04030	−0.0402958
	0.47105	0.4710498	−0.09121	−0.0912147	−0.03377	−0.0337734
	0.64765	0.6476541	0.01359	0.0135921	−0.02506	−0.0250638
	0.37427	0.3742730	0.14947	0.1494679	−0.00050	−0.0004986
	−0.28059	−0.2805924	−0.00354	−0.0035408	0.00006	0.0000633
	−0.10107	−0.1010694	−0.00573	−0.0057344	0.00011	0.0001124
	−0.01280	−0.0128006	−0.00223	−0.0022296	0.00002	0.0000172
	−0.00132	−0.0013249	−0.00036	−0.0003618

Table 9

Computational results relevant to phase velocities for a two-layer model.

T (s)	c (m∙s−1)		∂c/∂β1		∂c/∂β2
	Novotný	This paper	Novotný	This paper	Novotný	This paper
20	3441.81	3441.8133	0.74250	0.7425026	0.15013	0.1501333
30	3755.73	3755.7307	0.33374	0.3337405	0.47141	0.4714082
40	3896.68	3896.6754	0.14381	0.1438077	0.59832	0.5983238
			∂c/∂α1		∂c/∂α2
			Novotný	This paper	Novotný	This paper
20			0.14041	0.1404094	0.00231	0.0023095
30			0.16901	0.1690098	0.01454	0.0145360
40			0.13987	0.1398705	0.02492	0.0249162
			∂c/∂h		∂c/∂ρ
			Novotný	This paper	Novotný	This paper
20			−0.01981	−0.0198137	0.18560	0.1856024
30			−0.01897	−0.0189727	0.39702	0.3970193
40			−0.00965	−0.0096475	0.34035	0.3403527

Table 10

Computational results relevant to group velocities for a two-layer model.

T (s)	U (m∙s−1)		∂U/∂β1		∂U/∂β2
	Novotný	This paper	Novotný	This paper	Novotný	This paper
20	2864.63	2864.6298	1.09993	1.0999266	−0.28108	−0.2810790
30	3191.45	3191.4546	1.05173	1.0517331	−0.01968	−0.0196826
40	3585.94	3585.9400	0.45706	0.4570642	0.36729	0.3672883
			∂U/∂α 1		∂U/∂α 2
			Novotný	This paper	Novotný	This paper
20			0.02799	0.0279941	−0.00754	−0.0075417
30			0.21034	0.2103442	−0.01419	−0.0141852
40			0.23521	0.2352072	−0.00267	−0.0026709
			∂U/∂h		∂U/∂ρ
			Novotný	This paper	Novotný	This paper
20			0.00492	0.0049175	−0.27858	−0.2785814
30			−0.04428	−0.0442750	0.35718	0.3571744
40			−0.03018	−0.0301845	0.60146	0.6014577

Table 11

A six-layer model from Xia [1].

Layer	P-Wave Velocity (m∙s−1)	S-Wave Velocity (m∙s−1)	Density (kg∙m−3)	Thickness (m)
1	650	194	1820	2.0
2	750	270	1860	2.3
3	1400	367	1910	2.5
4	1800	485	1960	2.8
5	2150	603	2020	3.2
6	2800	740	2090	∞

Table 12

PDs of the fundamental-mode phase velocity with respect to β.

Frequency (Hz): [5, 10, 15, 20, 25, 30](Xia et al. 1999) [1]

$\left(\begin{array}{cccccc}0.018& 0.018& 0.022& 0.021& 0.017& 0.872\\ 0.130& 0.106& 0.062& 0.025& 0.022& 0.766\\ 1.067& 0.925& 0.313& 0.034& 0.017& 0.262\\ 0.155& 1.037& 0.967& 0.457& 0.145& 0.040\\ 0.293& 1.072& 0.517& 0.102& 0.012& 0.001\\ 0.520& 0.923& 0.202& 0.016& 0.000& 0.000\end{array}\right)$

(Lai & Rix 1998) [6]

$\left(\begin{array}{cccccc}0.01443& 0.01300& 0.01383& 0.01193& 0.01382& 0.86939\\ 0.12161& 0.08636& 0.04675& 0.02314& 0.01812& 0.76576\\ 0.95916& 0.84527& 0.28894& 0.02847& 0.01531& 0.25958\\ 0.14780& 0.96370& 0.92731& 0.45258& 0.13932& 0.03864\\ 0.28708& 1.04603& 0.49245& 0.10321& 0.01138& 0.00070\\ 0.51321& 0.89495& 0.20270& 0.01552& 0.00060& 0.00001\end{array}\right)$

(Cercato 2007) [35]

$\left(\begin{array}{cccccc}0.01809& 0.01834& 0.02219& 0.02036& 0.01750& 0.87242\\ 0.13002& 0.10646& 0.06174& 0.02467& 0.02225& 0.76580\\ 1.06766& 0.92490& 0.31304& 0.03359& 0.01665& 0.26204\\ 0.15460& 1.03665& 0.96729& 0.45739& 0.14507& 0.04024\\ 0.29284& 1.07203& 0.51690& 0.10261& 0.01137& 0.00074\\ 0.52024& 0.92354& 0.20196& 0.01592& 0.00060& 0.00001\end{array}\right)$

(Wu et al. this paper)

$\left(\begin{array}{cccccc}0.0180908& 0.0183407& 0.0221925& 0.0203622& 0.0174990& 0.8724199\\ 0.1300197& 0.1064578& 0.0617429& 0.0246696& 0.0222537& 0.7658011\\ 1.0676629& 0.9249004& 0.3130410& 0.0335878& 0.0166527& 0.2620402\\ 0.1546003& 1.0366473& 0.9672942& 0.4573936& 0.1450729& 0.0402373\\ 0.2928380& 1.0720290& 0.5168966& 0.1026072& 0.0113722& 0.0007421\\ 0.5202410& 0.9235443& 0.2019644& 0.0159170& 0.0005995& 0.0000107\end{array}\right)$