Surface waves are the most important attribute for discrimination between earthquakes and explosions. Traditionally, the M_s (Rayleigh wave magnitude) is measured at or around 20 sec period, to be used for (M_s: m_b) method of discrimination. The first estimate for M_s was introduced by Gutenberg (Gutenberg, 1945)that was applicable for the 20 sec period only and distances greater than 25^o. This method was improved by Vanek et al known as Prague formula(Vanek, 1962). They included the time period term in the empirical relation of Gutenberg to use the periods in the vicinity of 20 seconds. Basham reported that a great improvement in identification using Rayleigh waves was achieved if the measurements for surface wave magnitudes were restricted to the shorter period (8–14 second) and the epicentral paths were continental only(Basham, 1969a, b). Basham improved the Venek empirical relation for M_s for distances smaller than 25^o. That relation was further improved by Marshal and Basham for the period range of 8–14 sec(Marshall and Basham, 1972). Nuttli and Kim gave the empirical relation for the same conditions as in empirical relation of Marshal and Basham except that time period was around 20 sec. In all these works, empirical relations were developed and improved based on Rayleigh waves taken from the vertical component data. After the advent of broadband seismic systems and method of coordinate system rotation, it became possible to work on M_s from Love and Rayleigh waves separately. Surface wave magnitude widely used for yield estimation of explosions and discrimination between earthquakes and explosions in many studies like(Alewine III, 1972; Alexander and Rabenstine, 1967; Anglin, 1971; Anglin and Israelson, 1973; Basham, 1969a, b; Basham and Anglin, 1973; Bonner et al., 2008; Bonner et al., 2003; Bonner et al., 2006; Bonner et al., 2011; Chun et al., 2011; Gaber et al., 2017; Harkrider et al., 1994; McGarr, 1969; Murphy et al., 2013; Rezapour and Pearce, 1998; Russell, 2006; Selby, 2010; Shin et al., 2010; Stevens and McLaughlin, 2001; Zhang and Wen, 2013; Zhao et al., 2016; Zhao et al., 2008; Zhao et al., 2012; Zhao et al., 2014).

Seismic discrimination of distinguishing natural tectonic earthquakes from anthropogenic sources (such as underground nuclear explosions) remains foundational for seimic monitoring and verification regime. A classical and widely used discriminant is the comparison between a body wave magnitude (Mb) and a surface wave magnitude (Ms). The physical rationale is an isotropic explosion tends to radiate relatively more high frequency P-wave energy and less long-period surface wave energy. By contrast, a fault rupture earthquake typically couples more strongly into surface waves (Bonner et al., 2011; Nowroozi, 1986).

The Rezapour and Pearce used theoretical aspects of dispersion and geometrical spreading to develop a new surface wave magnitude equation.

Where A is the amplitude (in nanometers) of the instrument corrected ground motion at period of T second and Δ is the epicentral distance in degrees.

Bonner et al opposed Rezapour and Pearce for period less than 10 second and suggested that it needed path corrections. They found much better agreement between the Stevens and Murphy M_s values and the 7-sec modified Marshall and Basham estimates.Russell proposed a new surface-wave magnitude M_s(b) which differed from a traditional 20 second magnitude in that it used a phase matched Butterworth filter to measure a time-domain amplitude in a narrow band around any desired frequency, and then applied a correction for the source function. The main purpose of M_s(b) was to allow surface waves to be measured at regional distances at higher frequencies. The Russell’s formula is as follows

where is the amplitude of the Butterworth-filtered surface waves (zero-to-peak in nanometers) and is the filter frequency of a third order Butterworth bandpass filter with corner frequencies. This method can measure surface wave magnitudes at variable periods ranging from 8 to 25 second and at both regional and teleseismic distance.

Jessie L. Bonner and colleagues developed variable period, time domain surface wave magnitude measurement procedure applicable at both teleseismic and regional distances. In the context of discriminating surface waves (Love verses Rayleigh) using mb:Ms the Bonner technique is particularly relevant. By obtaining separate magnitudes for the Love and Rayleigh wave modes, one may compute ratios and thereby exploit differences in mode-excitation between earthquakes and explosions (Bonner et al 2010, 2011, 2013).

The complex crustal and tectonic environment of south Asia including the Pakistan region mb:Ms screening applied directly. Crustal heterogeneity, shallow source depths, path and site effects, anisotropy and regional attenuation all influence the surface wave mode content and amplitude. Therefore, the present study employs the Bonner technique to obtain separate Love and Rayleigh wave magnitudes, applies the refined mb:Ms, L and mb:Ms, R discriminants and evaluates them in the regional context. The study aims to (i) calibrate the discrimination boundary for regional data, (ii) assese whether Love vs Rayleigh ratio adds discrimination power beyond classic mb:Ms, and (iii) highlight implications for monitoring and verification in this tectonically active region (Bungum and Tjostheim, 1976; Krass, 2020).

In the present study, the broadband digital record of IRIS seismic station NIL (located at Nilore, Islamabad, Pakistan) has been used to check the differences in Rayleigh and Love wave magnitudes for the Earthquake events in and around Pakistan. A code has been written in Python language to estimate the Rayleigh and Love wave magnitudes at variable period range of 8–25 seconds. For this purpose, Russell’s magnitude formula (Eq. 02) is used. In this method, maximum amplitudes are found for variable periods using 3rd order Butterworth filter. For a single case, there are 18 readings of (maximum amplitudes in nanometers) Rayleigh waves and 18 readings of Love waves for periods ranging from 8 to 25 seconds. These maximum amplitudes and corresponding period will be used in Eq. 02 to find magnitude of Love or Rayleigh waves. The maximum magnitude among the 18 values will be taken as M_s(VMAX). For Rayleigh wave, vertical component data has been used while for Love wave the data of Transverse component has been used.

The data has been downloaded from IRIS Client FDSN using Obspy module of Python. Data has been decimated from its original sampling rate to 1 sample per sec, corrected for the instrument response, converted to displacement in nanometers, and rotated to transverse, radial, and vertical components. Rotated data has been used as input for our 2nd python code written for finding M_s(VMAX) for Love and Rayleigh waves. For earthquake study, NIL station data for the last ten years (2014 to 2023) has been used. Five regions have been investigated, which include Pakistan, India, Afghanistan, Iran and China. India region include events with depth less than or equal to 30 km and magnitude greater than or equal to 5for this study whilefor explosion studies, Chinese, North Korean, Indian and Pakistani events have been analyzed (Figure. 1). As for nuclear explosions only limited data was available on NIL station, so data from different international seismic stations with varying epicentral distances (800 to 5000 km) have been used to check our code for discrimination.

The analysis of broadband waveforms from both tectonic earthquakes and suspectedexplosion events reveals distinct characteristics in the Rayleigh and Love wave amplitudes.Tectonic earthquakes consistently display relatively stronger Love wave amplitudes,reflecting the shear-faulting mechanism, whereas explosion events produce dominantRayleigh waves due to isotropic source radiation. Eq. (2) was applied for the whole datasets of earthquakes occurred during a 10 years’ period from 2014 to 2023 (Figure. 1). The comparison between ofRayleigh and Love waves for all earthquake was estimated fortotal of 349 earthquakes from Afghanistan, China, India, Iran and Pakistan. Overall 89.4% results are above the line showing of Love waves are higher than or equal to that of Rayleigh waves (Table 1, Fig. 2). Individually for each region, the accuracyshown a reliable application of this method to discriminate earthquakes. The dominant time period for Rayleigh wave for all five regions lies in the range 11–19 sec and for Love wave 16–18 sec except Afghanistan region with low for both Rayleigh and Love waves (Fig. 3) (Nuttli and Kim, 1975). The relationship between the left- and right-window amplitude estimates was examined to evaluate the reliability of the VMAX method in discriminating seismic events based on surface wave characteristics.A consistent linear trend close to the 𝑀𝑏:𝑀𝑠 reference line was observed across all regions, indicating the stability of the VMAX-derived magnitudes. However, distinct regional differences were evident in the proportion of data points falling above or below the equality line.

A consistent linear trend close to the 𝑀𝑏:𝑀𝑠 reference line was observed across all regions, indicating the stability of the VMAX-derived magnitudes. However, distinct regional differences were evident in the proportion of data points falling above or below the equality line. In Afghanistan and Pakistan, 94.4% and 92.0% of the events, respectively, showed 𝑀𝑠(𝑉𝑀𝐴𝑋)𝐿 ≥ 𝑀𝑠(𝑉𝑀𝐴𝑋)𝑅, suggesting that the long-period energy recorded in the later time window was relatively stronger. This dominance of

𝑀𝑠(𝑉𝑀𝐴𝑋)𝐿 indicates higher energy release at longer periods, characteristic of natural tectonic earthquakes. Conversely, in India, Iran, and China, the majority of events (86.4%, 87.7%, and 86.4%, respectively) displayed 𝑀𝑠(𝑉𝑀𝐴𝑋)𝑅<𝑀𝑠(𝑉𝑀𝐴𝑋)𝐿, meaning that the right-window amplitude (earlier arrival) was lower than the later one. This trend corresponds to the nature of surface wave propagation in regions where structural heterogeneity and crustal scattering significantly affect the amplitude build-up in the later portions of the waveform. The strong clustering of earthquake data (red dots) above the reference line compared to the nuclear explosion data (black dots) indicates that natural earthquakes generate larger Rayleigh-wave amplitudes for the same body-wave magnitude, whereas explosions tend to produce relatively smaller 𝑀𝑠(𝑉𝑀𝐴𝑋) values due to their shallower and more isotropic source mechanisms.

Representative examples of surface wave recordings for an event from Pakistan (Fig. 4) further illustrate the time–frequency dependence of the VMAX-derived magnitudes. The event, with a body-wave magnitude 𝑀𝑏=4.8, was recorded at four stations (AAK, RAYN, MAKZ, and KURK) at epicentral distances ranging from ~ 1750 to ~ 2700 km. The recorded seismograms show a clear development of Rayleigh wave trains with varying dominant periods between 13–20 s.

At each station, the red and green waveforms correspond to amplitude measurements from two different period bands. The maximum amplitudes are observed in the range of 15–20 s, indicating that the surface wave energy is concentrated in the long-period range. Stations RAYN and KURK, situated at greater distances, display higher amplitude stability and clearer waveform envelopes, suggesting less attenuation and higher signal-to-noise ratio in the long-period bands. This stable amplitude pattern is reflected in the comparable 𝑀𝑠(𝑉𝑀𝐴𝑋) estimates across stations, demonstrating the consistency of the VMAX technique in identifying the maximum surface wave amplitude independent of moderate path effects.

Furthermore, the increase in amplitude with longer periods signifies that earthquake sources radiate more low-frequency energy compared to nuclear explosions, which typically exhibit stronger high-frequency components and shorter duration waveforms. The observed amplitude build-up in the 17–20 s band provides a strong indication of a tectonic source mechanism. These observations emphasize that the time–frequency representation of surface wave amplitudes can effectively differentiate between natural and artificial seismic sources.

Figure 5 summarizes the variation of 𝑀𝑠(𝑉𝑀𝐴𝑋) with period for selected events from North Korea, India, Pakistan, and China, as recorded at the TLY station. Each curve represents the magnitude variation as a function of time period (8–24 s), allowing direct comparison of how the VMAX amplitude changes across different source types and propagation paths. The North Korean nuclear explosion (1993 event, 𝑀𝑏=5.9 shows a systematically higher 𝑀𝑠(𝑉𝑀𝐴𝑋)𝑅 value (4.29) compared to 𝑀𝑠(𝑉𝑀𝐴𝑋)𝐿 (3.77), yet both magnitudes are lower than those for natural earthquakes of similar size. This reduction in 𝑀𝑠is consistent with the known behavior of explosion-generated surface waves, which are relatively weak due to the shallow, isotropic source and lack of significant shear displacement. The relatively flat curve of 𝑀𝑠(𝑉𝑀𝐴𝑋) across increasing periods also indicates that explosion energy does not preferentially radiate into longer periods, confirming its distinct spectral signature. For the Indian and Pakistani earthquakes, the 𝑀𝑠(𝑉𝑀𝐴𝑋) curves exhibit clear period-dependent variations. In both cases, 𝑀𝑠(𝑉𝑀𝐴𝑋) reaches a maximum in the range of 15–18 s and then slightly decreases toward 24 s. This gradual rise and fall correspond to the natural period of Rayleigh wave generation in regional crustal structures. The Pakistani event (𝑀𝑏=4.8) shows 𝑀𝑠(𝑉𝑀𝐴𝑋)𝑅=3.38 and 𝑀𝑠(𝑉𝑀𝐴𝑋)𝐿=3.22, indicating that both left- and right-window amplitudes are nearly equivalent, suggesting a balanced spectral energy distribution typical of an earthquake source. The Chinese event (𝑀𝑏=5.1) shows 𝑀𝑠(𝑉𝑀𝐴𝑋)𝑅=4.14 and 𝑀𝑠(𝑉𝑀𝐴𝑋)𝐿=3.82, which are higher than those for Pakistan and India, reflecting the effect of regional propagation path and possible crustal amplification near the TLY station. The consistent amplitude growth with increasing period further supports the earthquake source origin, where shear dislocation contributes strongly to surface wave generation.

The comparative behavior of 𝑀𝑠(𝑉𝑀𝐴𝑋) across regions and event types confirms that the VMAX approach provides a robust metric for distinguishing earthquakes from underground nuclear explosions. Explosion events are characterized by (1) systematically lower 𝑀𝑠(𝑉𝑀𝐴𝑋) values, (2) flatter or decreasing trends with period, and (3) smaller amplitude differences between left and right time windows. Earthquakes, in contrast, show (1) stronger period dependence, (2) larger long-period amplitudes, and (3) a tendency for 𝑀𝑠(𝑉𝑀𝐴𝑋)𝐿 to exceed 𝑀𝑠(𝑉𝑀𝐴𝑋)𝑅. The high percentage of 𝑀𝑠(𝑉𝑀𝐴𝑋)𝐿>𝑀𝑠(𝑉𝑀𝐴𝑋)𝑅 in Afghanistan and Pakistan (over 90%) reflects that regional earthquakes consistently radiate significant long-period surface wave energy, whereas the reduced percentages in India, Iran, and China may be attributed to crustal heterogeneity, attenuation, or site effects influencing the amplitude balance. The combined waveform and magnitude analyses therefore demonstrate that the VMAX method can effectively capture subtle differences in surface wave characteristics that correlate with source type and propagation path.These findings validate the effectiveness of the VMAX approach as a practical and stable tool for regional seismic discrimination and event characterization.

This study confirms the reliability of the VMAX-based Love–Rayleigh wave discriminant for distinguishing earthquakes from underground nuclear explosions in and around the Pakistan region. Using broadband NIL station data (2014–2023), the method correctly identified 89.4% of 349 shallow earthquakes, with regional accuracies up to 94.4%, and effectively captured period dependence across 8–25 s, highlighting the advantage of variable-period analysis.

For explosion events from North Korea, China, India, and Pakistan, the method achieved 88.3% overall accuracy, performing best at regional distances (< 2000 km) where surface-wave amplitudes are strong and coherent. Larger events yielded clearer discrimination, while small-magnitude tests were more affected by noise. Overall, the VMAX discriminant provides a simple, data-efficient, and regionally adaptable approach for seismic source characterization. Its integration into regional monitoring networks can enhance nuclear test verification and earthquake explosion discrimination. Future improvements should focus on attenuation corrections and multi-station analyses to refine performance at teleseismic distances.