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

: 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 ﬁrst-order partial di ﬀ erentiations and approximates the mixed second-order partial di ﬀ erentiation term with a central di ﬀ erence scheme. We provide examples to demonstrate the e ﬀ ectiveness 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.


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 phaseand/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).

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 α m , the S-wave velocity β m , the density ρ 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: where ω is the angular frequency, c is the phase velocity, m is the sequence number for layers, and n + 1 indicates the half-space. X * m is a row vector with different lengths depending on which wave (Love or Rayleigh wave) is concerned. T * m and 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.
Appl. Sci. 2019, 9,5214 3 of 21 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 surfacewave 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 hemianalytical 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).

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α , the S-wave velocity , the density and the thicknessℎ . 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,2, ⋯ , where ω is the angular frequency, c is the phase velocity, m is the sequence number for layers, and n+1 indicates the half-space. * is a row vector with different lengths depending on which wave (Love or Rayleigh wave) is concerned. * and * 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.

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

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, for the case of Rayleigh waves, or for the case of Love waves. Hence, the completed form of F is Suppose we want to compute ∂c/∂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, ∂c/∂p t can be written as: In the following, let us first consider ∂F/∂c. On the basis of recursion (1) and the derivative rules, the expression for ∂F/∂c can be derived as follows: Obviously, the calculation of ∂V * n+1 /∂c is straightforward, whereas ∂X * n+1 /∂c has to be calculated recursively starting from the first layer via formula (1): For a clearer description of the algorithm (Table 1), let us introduce a new vector, i.e., Y * m = ∂X * m /∂c, thus (6) can be reformulated as follows: 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 ∂V * n+1 /∂c, ∂T * m /∂c, ∂V * n+1 /∂p n+1 , and ∂T * t /∂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. Table 2. Pseudo-code for computation of ∂F/∂p t .

Pseudo-code 2
The computation of ∂F/∂p t Input: angular frequency ω, phase velocity c, earth model Q R or Q L , and the sequence number t of the objective layer 1: Initialize one vector and one matrix: for m = 1-n do 10: if m equals to t then 11: Calculate

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]: 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.

Fully Analytical Method
The first order derivative ∂c/∂ω in (8) also can be calculated by using the implicit function theory: Although the pseudo-code in Table 1 is for ∂F/∂c, it is apparently applicable to ∂F/∂ω, 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]: The second-order mixed partial derivative in (10) can be derived by further differentiating (4) with respect to ω [31]: In the following, we take ∂ 2 F ∂c 2 as an example to illustrate how to calculate the second-order partial derivatives analytically. As for ∂ω∂c , they can be calculated in the same way. The computational procedure for ∂F/∂c (Table 1) can be performed in a recursive manner: Therefore, we differentiate the last expression above with respect to c and thus obtain Then, in order to obtain ∂Y * n+1 /∂c, we shall resort to the second expression in (12): Similarly, if we introduce a new vector Z * m = ∂Y * m /∂c, then (14) can be reformulated as follows: Obviously, 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 ∂ 2 F/∂c 2 (Table 3). See Appendix C for the elements of ∂ 2 T * m /∂c 2 and of ∂ 2 V * n+1 /∂c 2 for Love waves. Table 3. Pseudo-code for computation of ∂ 2 F/∂c 2 . update xold: xold = xnew 8:

Hemi-Analytical Method
In the fully analytical method, the crucial second-order partial derivative term ∂ 2 c ∂ω∂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 ∂U/∂p t is desired at the angular frequency ω 0 and that the phase velocity c 0 and group velocity U 0 are computed at ω 0 . Define two other frequencies around ω 0 : ω +1 = ω 0 e +δ and ω −1 = ω 0 e −δ . δ is a small quantity, denoting a perturbation of log ω; see [36] for a more detailed explanation. And suppose the quantities c +1 , U +1 , and ∂c +1 /∂p t are corresponding to ω +1 and have been obtained. The same hypothesis is applied to the quantities c −1 , U −1 , and ∂c −1 /∂p t at ω −1 . Then the PD of the group velocity at ω 0 can be computed via 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 ∂ 2 c ∂ω∂p t approximated by the central finite-difference.

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 T * m and of dispersion function F. To make it clearer, we list these files in Tables A1 and 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 ∂F/∂β 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.

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 × ∂c/∂h and 10 × ∂U/∂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 ρ = ρ 2 /ρ 1 and plotted the curves of ∂c/∂ρ and ∂U/∂ρ. Therefore, to obtain the same curves, we have made use of the chain rule:  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.    Obviously, the computational results in Figure 2 are rather consistent with the work of Novotný [31] (see Figures 1-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 Tables 1 and 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 Tables 7 and 8 for detailed comparisons. It should be pointed out that Novotný [30] (Figures 1-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.

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 (α 1 = 6000 m·s −1 and α 2 = 8000 m·s −1 ) [32]. Figure 4 displays the related computational outcomes of fundamental-mode phase and group velocity. They show apparently good agreements with [32] (Figures 1-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.

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] (Figures 3-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.      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 10.

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.
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)"

Appendix B
The expressions listed below are used for the construction of the first-order PDs of matrices T * m and V * n+1 for Love waves. For clarity, we use the same notation as Cercato [35], i.e., the prime mark (e.g., s ) 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 β = ∂s/∂β and C β β = ∂C β /∂β).