Efficient Computation of Dispersion Curves in Low‐Velocity‐Layered Half‐Spaces

The computation of dispersion curves for surface waves holds significant importance across various applications, such as the simulation of displacement in surface waves (Harkrider, 1964; Kerry, 1981; Zeng and Anderson, 1995). In addition, it plays a crucial role in the interpretation of Earth structures, as evidenced by studies conducted by Xia et al. (2003), Shapiro et al. (2005), Luo et al. (2007), and Pan et al. (2022).

The methods employed for determining surface‐wave dispersion curves in a flat‐layered Earth model originated from the work of Thomson (1950) and Haskell (1953) in the early 1950s. Initially, they utilized a matrix approach to address the eigenvalue problem inherent in the differential equations system. However, the original Thomson–Haskell formalism proved to be unstable and encountered various numerical challenges, such as issues with numerical overflow and precision loss at higher frequencies (Buchen and Ben‐Hador, 1996). Subsequent enhancements to the formalism were introduced to mitigate these issues, including the delta matrix (Pestel et al., 1964; Dunkin, 1965), the reduced delta matrix proposed by Watson (1970), the Schwab–Knopoff method by Schwab and Knopoff (1970, 1972), the Abo‐Zena method by Abo‐Zena (1979), and the fast delta matrix technique developed by Buchen and Ben‐Hador (1996).

Differing from the Haskell matrix method, the method of the generalized reflection and transmission (GRT) was first presented by Kennett (1983) and Luco and Apsel (1983), subsequently improved by Chen (1993) and Hisada (1994, 1995). The iterations by Chen and Hisada have demonstrated the capability to yield stable and precise phase velocities across a wide range of frequencies, encompassing both low and high values. Pei et al. (2008) introduced a rapid GRT approach aimed at enhancing computational efficiency through the direct computation of GRT coefficients. Building upon this, Wu and Chen (2016, 2022) developed a flexible normal modes solver utilizing the GRT method specifically designed for intricate horizontal stratified models, which may incorporate complex features such as low‐velocity layers and fluid interlayers.

A variety of approaches have been advanced for determining dispersion curves, with indications of the prevalent presence of the mode‐missing issue, particularly in models incorporating low‐velocity layers. Research has demonstrated the existence of a phenomenon known as “mode‐kissing” in layered media with low‐velocity layers, for which two distinct dispersion curves exhibit differing vibrational energy at kissing points (Liu et al., 2009; Gao et al., 2016). The instances of “mode‐kissing” typically coincide with locations where missing roots are identified. He et al. (2006) introduced the notion of a set of secular functions to address the computational challenge of identifying the trapped modes associated with energy concentrated predominantly in the low‐velocity layer. Furthermore, the identification of missing roots can be achieved through the comparison of roots at adjacent frequencies, as demonstrated by Wu and Chen (2022), to mitigate the issue of missing roots. However, the efficiency challenge of the GRT method remains a significant hurdle (Buchen and Ben‐Hador, 1996), particularly in instances for which a considerable number of forwardings are required for surface‐wave inversion.

In the field of surface‐wave inversion, two predominant methodologies for inversion are frequently employed. The first approach involves the pursuit of a set of unknown parameters through the optimization of a data misfit function (Sambridge, 1999; Pan et al., 2022). The second approach, which is increasingly favored, entails the probabilistic sampling of the posterior probability density function within a Bayesian framework (Bodin et al., 2012; Shen et al., 2013; Sambridge, 2014). This Bayesian approach offers significant advantages over the optimization‐based method, primarily due to its reduced susceptibility to convergence within local minima.

However, it is important to acknowledge the trade‐offs associated with the Bayesian methodology. Specifically, the accurate estimation of the posterior probability density often necessitates a considerable number of evaluations of dispersion curves, which can reach into the millions throughout the iterations of the Markov chain. Furthermore, the models employed for the computation of dispersion curves are typically generated randomly, which may result in complex representations of models. Consequently, it becomes imperative that the forward modeling algorithm utilized is both stable across a wide range of scenarios and efficient enough to accommodate the demands of extensive evaluations.

In light of the identified application requirements and obstacles, a novel approach is introduced for computing dispersion curves by applying an adaptive root‐finding procedure. The algorithm is designed to introduce supplementary samples between adjacent modes in the vicinity of kissing points. This is achieved by sequentially calculating the roots for a degraded model, which involves the systematic removal of the low‐velocity layer in conjunction with the underlying layers. Following the implementation of the root‐searching scheme, comprehensive dispersion curves can be derived.

The computation of dispersion curves primarily comprises two fundamental components: the evaluation of the dispersion function and the specification of sampling points on phase velocities for assessment. The first component involves defining the dispersion function, which serves to ascertain whether a given phase velocity corresponds to the dispersion curves. The algorithm utilized in this study is the delta matrix method, which is also employed in the Computer Programs in Seismology (CPS; Herrmann, 2013), thereby facilitating a comparative analysis. It is noteworthy that the dispersion function at a single frequency exhibits multiple roots, with varying distances between adjacent roots. Because of the multiroot characteristics of the dispersion function, it is imperative to identify appropriate locations for evaluating the dispersion function, subsequently determining the phase velocities that satisfy the condition that the dispersion function equals zero. The primary focus of this article is on the latter aspect, whereas the selection of the dispersion function itself lies outside the scope of this investigation.

The delta matrix method (Dunkin, 1965) is employed in this study to evaluate the dispersion function.

To elucidate the distribution pattern of the dispersion function, a four‐layer model incorporating a low‐velocity layer was utilized. The parameters of the model are detailed in Table 1. The dispersion function D(c,f) depicted as black curves in Figure 1 is scrutinized across varying phase velocities at both 10 and 30 Hz frequencies. The dispersion function has undergone a sign operation applied to its values to enhance the clarity of its presentation. Phase velocities that correspond to solutions of the dispersion equation (equation 1) align with the theoretical dispersion curves of interest, which can be efficiently determined through the application of the bisection algorithm (Burden and Faires, 2016), a root‐finding method that iteratively narrows bracketing intervals by halving.

At a specific frequency, the dispersion equation demonstrates multiple roots that correspond to the phase velocities of dispersion curves. Dispersion is commonly ascertained by exploring a range of phase velocities that span the potential extremes, which are contingent upon the lowest shear‐wave velocity within the layers and the shear‐wave velocity of the half‐space. A straightforward approach involves conducting a search using a constant increment $δc$⁠. However, to mitigate the phenomenon of mode jumping, it is necessary to utilize a significantly small value for $δc$⁠, which consequently leads to an increase in computational time due to the requisite numerous evaluations of the dispersion function. The sdisp96 program from CPS addresses this issue by implementing a small $δc$ at high frequencies and subsequently endeavors to track the dispersion of a mode as it transitions to lower frequencies. Furthermore, sdisp96 employs a more refined search strategy within the region where mode jumping is likely to occur. The underlying code for this process is relatively complicated. In contrast, the alternative approach proposed in this article involves the application of an efficient procedure at each frequency, leveraging recent findings (Fan et al., 2007; Ke et al., 2011; Jiang et al., 2024) on the relationship between dispersion curves and velocity models to minimize the number of functional evaluations required.

The anticipated modal quantities determined by $N(c)$ for model 1 at 10, 50, and 90 Hz frequencies have been calculated and are depicted as black lines in Figure 2a, demonstrating a strong agreement with the enumeration of modes from the theoretical dispersion curves. Through the process of determining the roots of equation (4) and subsequent refinement, samples of the phase velocity used in root finding at three frequencies can be acquired and are illustrated in Figure 2b–d, respectively. The parameter $ϵ$ is assigned a value of 0.5, whereas M is designated as 2 in this article.

The establishment of predetermined phase‐velocity samplings for subsequent root‐finding processes is articulated as follows.

1. Utilizing the established velocity model, it is essential to delineate the maximum and minimum potential values of phase velocity associated with the specific wave type. Typically, the maximum phase velocity, denoted as $cmax$⁠, is represented by the shear‐wave velocity in the half‐space $βhalf−space$⁠. Conversely, the minimum phase velocity $cmin$ is conventionally determined as the smallest shear‐wave velocity within the model, with the exception of Rayleigh waves, for which it corresponds to the Rayleigh‐wave phase velocity linked to $β1$⁠.

2. An estimation of the maximum number of modes $N(cmax)$ can be derived from equation (2).

3. Subsequently, a search should be conducted over the range $c∈[cmin,cmax]$ to identify $cj$ utilizing equation (4). The value $N(cmax)$ is instrumental in establishing an initial phase velocity increment for this search. To enhance the list of potential phase velocities, apply the method of interval halving between values a specified number of times, denoted as M.

4. In addition, the list may be further expanded by incorporating a fraction of the minimum shear‐wave velocity $βmin$ for instance, by including $0.9βmin$⁠. In the context of Rayleigh waves, it is pertinent to also include the Rayleigh‐wave phase velocity associated with $β1$⁠.

The predetermined sampling protocol is implemented in the $get_samples$ function of the $src/disp.cc$ source file within the mode‐kissing codebase (see Data and Resources).

The intervals housing the roots of the dispersion function are determined by evaluating D(c, f) at adjacent initial sampling points, which are then employed for root localization through the bisection algorithm. The pseudocode detailing the root‐finding process based on the initial phase‐velocity samplings in Algorithm 1. The root‐finding process resides in the $search$ function of the same source file (⁠$src/disp.cc$⁠).

Utilizing the root‐searching scheme, dispersion curves for model 1 can be computed using the fast delta matrix algorithm. The dispersion curves are determined up to a frequency of 100 Hz with an interval of 0.1 Hz, and up to a maximum of 31 modes. As illustrated in Figure 2, multiple absent roots are situated at locations where “mode‐kissing” phenomena occur. To explain what happens at “mode‐kissing” points, the dispersion function at 19.7 Hz is studied, which is the location of the blue region in Figure 3. Figure 4 illustrates the dispersion function. The red‐vertical lines delineate the phase velocities corresponding to theoretical dispersion curves, whereas the black‐vertical lines signify the samplings determined through the root‐searching algorithm. Notably, the absence of roots is observed at a phase velocity approximately equal to 0.267 km/s. Upon closer examination of this region, it is evident that there are no samples positioned between the two roots, consequently resulting in the absence of two roots that should have existed.

Liu et al. (2009) discussed the phenomenon of “crossing” between two distinct modes of dispersion curves in layered media containing a low‐velocity layer. They designated the mode in which the vibrational energy is concentrated primarily near the surface of the medium as the R mode and referring to the mode for which the energy is confined within the low‐velocity layer as the $S2$ mode. In the vicinity of the cross points, a transformation between the R mode and the $S2$ mode occurs through coupled modes along each curve. Gao et al. (2016) conducted an analysis on the phenomenon of “mode‐kissing” by examining the eigendisplacements’ characteristics. They observed that the Rayleigh‐wave mode in proximity to the kissing points undergoes a transition in its type, leading to a scenario for which a single Rayleigh‐wave mode encompasses multiple mode types. This conversion of mode types results in the manifestation of the “mode‐kissing” phenomenon in dispersion images.

The eigendisplacements for model 1 are calculated using the GRT coefficient method as described by Chen (1993). The eigendisplacements corresponding to frequencies within the blue region depicted in Figure 3 are visualized in Figure 5. Specifically, Figure 5a illustrates the eigendisplacements associated with the first mode of dispersion curves, whereas Figure 5b depicts those for the second mode. Through an analysis based on eigendisplacements, it can be observed that the energy distribution of the right segment of the “crossed” frequency pertaining to the first mode and the left segment corresponding to the second mode is predominantly concentrated within the first layer. Conversely, the energy distribution of the left segment of the “crossed” frequency associated with the first mode and the right segment linked to the second mode is primarily concentrated within the low‐velocity layer. The concept of “mode” involves the enumeration of phase velocities ranging from lower to higher values at a specific frequency. However, it is noteworthy that an alternate definition of mode can be conceptualized from an energy standpoint. In practical scenarios, for which receivers are typically positioned on the surface, dispersion curves indicative of energy concentration in close proximity to the surface are readily discernible (Gao et al., 2016).

As a consequence of the case, He et al. (2006) introduced a methodology grounded in the GRT coefficient method utilizing a family of secular functions, as opposed to a single secular function, to address the issue of missing roots in the computation of dispersion curves for models incorporating low‐velocity layers. The secular function is conceptually analogous to the dispersion function delineated in this study. Each secular function within a set is specifically delineated for distinct layers. Inspired by this concept, we propose a novel algorithm for the computation of dispersion curves with additional samplings targeted at “mode‐kissing” points.

To address the challenge of missing roots, various solutions have been suggested. Among these, a direct approach involves enhancing the density of samples such that there is a sampling point positioned between every two adjacent roots at the intersection point. However, a limitation of this method arises due to the minimal distance between two roots, particularly evident at higher frequencies, necessitating a substantial number of samples to be taken to mitigate the occurrence of missing roots. An alternative method involves verifying the presence of absent roots at neighboring frequencies (Wu and Chen, 2022). This approach necessitates a relatively reduced frequency spacing to ensure timely detection of any potentially missing modes prior to solver resolution. Herrmann (2013) employed a jumping technique, which involves detecting the relationship of roots between a specific frequency and 99% of its value. It is imperative, from an intuitive standpoint, to generate a sample at the intersection point between two adjacent roots because this represents a primary objective of our algorithm.

We formulate a model that excludes the low‐velocity layer and the subsequent layers; specifically, for model 1, this entails incorporating only the first two layers. The dispersion curves corresponding to these models are presented in Figure 6. Notably, the dispersion curves of the degraded model intersect the mode‐kissing points of the original dispersion curves, a finding that is consistent with the eigenfunction analysis of the mode‐kissing phenomenon discussed in the preceding section.

Furthermore, we provide an enhanced view of specific locations where missing roots occur, as illustrated in Figure 3a. It is evident from Figure 7 that the dispersion points identified for the degraded model are situated between two modes of the dispersion curves of the complete model. This observation suggests that the dispersion points from the degraded model may serve as supplementary samples for identifying the roots of the complete model.

The proposed algorithm (Algorithm 2) computes dispersion curves through hierarchical model degradation. First, all low‐velocity layers are identified to construct the degraded model set $Mlvl$⁠, in which each $m∈Mlvl$ removes a shallow low‐velocity layer and its underlying layers. For every degraded model m:

1. initial phase velocity samples $initSamples$ are generated from equations (2) and (4),

2. these samples are augmented with roots addSamples discovered in previous iterations,

3. a root‐finding process (Algorithm 1) is applied to the combined sample set initSamples∪addSamples, and

4. newly found roots cList are appended to addSamples for subsequent iterations.

After processing all degraded models, the original undegraded model is analyzed: (1) final initSamples are generated with expanded velocity bounds, (2) all accumulated samples (initSamples∪addSamples) drive the root search, and (3) the complete solution cList is returned. This progressive refinement prevents mode omissions while maintaining computational efficiency through sample reuse across model variations. The high‐level workflow of multimodal root finding via degraded model referencing is orchestrated by the driver script $src/main_forward.cc$⁠.

The dispersion curves derived from the proposed algorithm are presented in Figure 8. These curves are computed up to a frequency of 100 Hz with an interval of 0.1 Hz. The original dispersion curves, obtained without supplementary sampling, are depicted in Figure 3a. The complexity associated with the computation of dispersion curves primarily stems from the frequent occurrence of root‐missing phenomena at intersection points, commonly referred to as mode‐kissing points. To address this issue, the algorithm delineated in Algorithm 2 employs strategic sampling between adjacent roots at these kissing points, thereby reducing the incidence of root missing.

To evaluate the efficacy of the computed dispersion curves, the results of the proposed methodology were compared with those generated by the widely used software application, CPS, as described by Herrmann (2013). In CPS, the parameter FACR (a predefined variable) within the sprep96 is utilized to control the root search process. Lower values of FACR facilitate faster computations but may result in the omission of higher modes. For the present study, the dispersion curves obtained using CPS were generated with the default FACR value of 5.0. Although increasing FACR can diminish the number of missing roots, it still yields several root missing even setting FACR to 20.0 and significantly prolong computation time.

The current study employs the proposed methodology to generate comprehensive dispersion curves, specifically addressing the absence of roots at high frequencies observed in the CPS program, as shown in Figure 8a. The proposed method successfully identifies a total of 16,042 dispersion points, in contrast to the 15,661 points yielded by the CPS program. In terms of computational efficiency, the proposed algorithm completes the calculation in 0.226 s, whereas the sdisp96 function in CPS requires 2.328 s, demonstrating approximately a tenfold increase in speed. Given that both the proposed method and CPS utilize the delta matrix method for computing the dispersion function, the primary factor influencing computation time is the selection of phase velocity sampling points for evaluating the dispersion function, which remains the focal point of improvement in this article.

Furthermore, we compare our algorithm with the program developed by Wu and Chen (2022), referred to as Surface Modes. The evaluation of its dispersion function is predicated upon the GRT method (Chen, 1993). The initial phase velocity samplings of SurfaceModes are generated in accordance with the algorithm established by Fan et al. (2007). To address the challenge of missing roots, the authors implement a verification process by comparing roots across adjacent frequencies. In instances where modes are missing, they identify the potential interval within which the lost root may exist and subsequently employ a root‐finding method again to locate it.

Utilizing the default parameters SurfaceModes provided, we present a comparison of the dispersion curves computed by our algorithm and SurfaceModes in Figure 8b. Notably, it is observed that SurfaceModes encounters a missing root. The computational time recorded for SurfaceModes is approximately 9.40 s after we have commented codes that output messages to stdout. The relatively prolonged execution time can be attributed to two primary factors. First, the GRT method used for evaluating the dispersion function is intrinsically slower compared to alternative methods, as indicated by Buchen and Ben‐Hador (1996). Second, to reduce the occurrences of missing roots, the authors have opted to increase the density of initial phase velocity samplings.

The numerical implementations of CPS and SurfaceModes employ Fortran, whereas the algorithm proposed in this study is developed in C++. Variations in computational time across these frameworks may arise, in part, from disparities in compiler optimization efficiencies between programming languages. However, this study focuses on advancing phase velocity sampling strategies for dispersion curve analysis, rather than directly comparing algorithmic implementations of dispersion function computations. To isolate the efficiency gains attributable to sampling methodology, we quantify computational demand using the total number of dispersion function evaluations—a metric independent of programming language or hardware.

For model 1, when resolving dispersion curves up to 100 Hz at 0.1 Hz increments, the proposed algorithm requires 425,075 dispersion function evaluations. By comparison, CPS necessitates 6,330,310 evaluations, and SurfaceModes requires 2,654,151 evaluations. These results demonstrate a reduction in computational load by a factor of 14.9 relative to CPS and 6.2 relative to SurfaceModes, confirming the superior efficiency of our phase velocity sampling approach.

In this study, we introduce a stable and effective approach for calculating the full dispersion curves of a multilayered model incorporating low‐velocity layers. The primary focus of our discourse pertains to models that incorporate low‐velocity layers; however, it is important to highlight that the methodology proposed herein is equally applicable to models that do not include such layers. The essence of our algorithm lies in its adeptness at generating samplings between adjacent roots, a phenomenon commonly referred to as “mode‐kissing.” Leveraging a root‐searching scheme and an algorithm designed for predicting cross points, our proposed method demonstrates superior efficiency for calculating dispersion curves of surface waves.

The computational architecture developed in this study directly enhances the feasibility of Bayesian surface‐wave inversion approaches that require massive evaluations of dispersion curves. Although Markov chain Monte Carlo methods provide superior uncertainty quantification compared to deterministic optimization (Bodin et al., 2012), their practical implementation has been historically constrained by two interdependent factors: (1) the computational burden of repeated forward modeling evaluations, and (2) the algorithmic stability across geologically plausible parameter spaces. Our methodology addresses this dual challenge through the numerical innovations, achieving an approximately 90% reduction in computation time relative to the CPS package (Herrmann, 2013) while maintaining robust results across the 0.1–100 Hz frequency spectrum.

This enhanced performance profile enables more rigorous exploration of posterior probability distributions, particularly critical in resolving nonunique solutions inherent to surface‐wave inversion problems. The improved forward modeling efficiency permits Bayesian sampling at higher‐dimensional parameterizations without prohibitive computational cost. Field data applications can now realistically incorporate complex prior distributions and transdimensional model parameterizations while adhering to practical runtime constraints. Future extensions could integrate machine‐learning‐accelerated likelihood evaluations to further push the boundaries of probabilistic inversion frameworks.

The model data utilized in this study have been sourced from the references cited in accordance within the article. All of the model data used in this study are from the references cited appropriately in the article. The data visualization is created through the utilization of the Matplotlib library (https://matplotlib.org/). The codes for the proposed method are available at https://github.com/pan3rock/mode-kissing. The codes for the Computer Programs in Seismology (CPS) are available at https://rbherrmann.github.io/ComputerProgramsSeismology/index.html. The codes for SurfaceModes are available at https://github.com/BranchInVine/SurfaceModes. All websites were last accessed in March 2025.

The authors acknowledge that there are no conflicts of interest recorded.

The authors extend their gratitude to the anonymous reviewers and Associate Editor Cezar I. Trifu for their rigorous critique and invaluable suggestions, which greatly enhanced the technical clarity and practical relevance of this work. This work is supported by Shenzhen Key Laboratory of Deep Offshore Oil and Gas Exploration Technology (Grant Number ZDSYS20190902093007855), Guangdong Provincial Key Laboratory of Geophysical High‐resolution Imaging Technology (Grant Number 2022B1212010002), and Shenzhen Science and Technology Program (Grant Number KQTD20170810111725321).