Theoretical Issues with Rayleigh Surface Waves and Geoelectrical Method Used for the Inversion of Near Surface Geophysical Structure

We numerically simulate the field measurements of Rayleigh surface waves and electrical resistivity in which the target depth is set to be less than 50-m. The Rayleigh surface waves are simulated in terms of fundamental mode group and phase velocities. The seismic field data is assumed to be collected through a conventional shot-gather. The group velocities are found from the application of the multiple filter technique in a single-station fashion while for the phase velocities the slant stacking, or linear radon transform are applied in fashion of multichannel analysis of surface waves (MASW). The average seismic structure from the source to the receiver (or geophone) is represented by the group velocity curve while the average seismic structure underneath the geophone array is represented by the phase velocity curve. The single-station group velocity curves are transformed into local group velocity curves by setting a linear system through grid points. The shear-wave velocity cross section underneath the examined area is constructed by inverting these local group velocity curves. The electrical resistivity structure of the underground is similarly studied. The field compilation of the resistivity data is assumed to be completed by the application of the multiple electrode PolePole array. The actual resistivity assemble underneath the analyzed area is inverted by considering the apparent (measured) resistivity values. Unique forms such as ore body, cavity, sinkhole, melt, salt, and fluid within the Earth may be examined by joint interpretation of electrical resistivities and seismic velocities. These formations may be better outlined by following their distinct signs such as high/low resistivities and high/low seismic velocities.


Introduction
The elastic P (compressional) and S (shear) waves created by either passive (e.g., an earthquake) or active (e.g., a sledgehammer) sources are regularly employed to study the seismic velocity structure within the Earth. The velocity structure, which we consider in near surface, is composed of low frequency component (i.e., multilayered composition with sub-horizontal layer interfaces) plus high frequency component made of mostly inclusions such as ore body, salt, melt, fluid, sinkhole, and cavity. Herein, we are principally interested in studying the low frequency multilayered component near the surface (<50 m). The depth range shallower than 30-m is particularly important to investigate the foundation properties characterized by 30 (shear-wave velocity) relevant to the earthquake mitigation studies [1][2][3]. Besides the body waves (i.e., P and S waves) there exist elastic surface waves, i.e., Rayleigh and Love surface waves.

184
The Rayleigh surface waves, which we currently consider, are created by constructive interference of multiply refracted and reflected compressional (P) and vertically polarized shear waves (SV). Because of the multilayered nature of the seismic velocity structure within the Earth, this constructive interference, which also includes wave-type conversions from SV-to-P and P-to-SV, produces dispersive waveforms recorded by seismic stations [4]. The seismic velocities normally increase with depth and there are cases of which some low velocity zones indicative of mechanically weak zones and partial melting develop in the stratified structure. On the other hand, emplacements of magmatic origin and ore body developments may cause high velocity zones.
The integrated geophysical interpretation of the sub-surface [5][6][7] is more advantageous than that one of these geophysical methods (i.e., seismic, electrical, gravity, magnetics, and electromagnetics) is applied alone. There is not a single geophysical method that is sensitive to all physical properties (e.g., velocity, attenuation, anisotropy, magnetic permeability, electrical conductivity, density) of the Earth. For instance, the surface waves are primarily sensitive to Swave (shear) velocities while the corresponding sensitivity to the density is only secondary. The gravity method combined with the seismic method can be used better to determine the density distribution within the Earth. Therefore, a combination of different geophysical methods may be useful to overcome the sensitivity problem, i.e., one geophysical method can help resolve ambiguities in another geophysical method [8][9]. In this respect, we consider the geoelectrical method alongside the Rayleigh surface waves in a joint manner to study the near surface Earth assembly. The geoelectrical method based on examining the electrical conductivity anomalies developing from variations from one rock unit to another one has found many uses in geophysics [10][11][12][13][14][15].
The surface waves (i.e., Rayleigh surface waves) have found many applications in studying the subsurface shearwave velocities underneath different locations [16][17][18][19][20][21]. The low frequency sub-surface is characterized by layer thickness (ℎ), shear-wave velocity ( ), compressional-wave velocity ( ), density ( ), shear-wave quality factor ( ), and compressional-wave quality factor ( ). We set Poisson's ratio to = 0.25 from which the is computed via the ⁄ ratio. In addition, we consider the electrical resistivity ( ) in each layer. Figure 1 illustrates the stratified one-dimensional (theoretical) model structure that we currently utilize. We compute the synthetic seismogram of Rayleigh surface wave to characterize the wave propagation media. The electrical current flow (direct current -DC) within the model structure is considered by computing the apparent resistivity pseudo-section. The DC surveys use current generator, voltmeter, and electrical contact (electrode) with the ground. The inversion algorithm is designed to recover the subsurface shear-wave velocities ( ) and electrical resistivities ( ) in each layer with a certain thickness (ℎ).

Figure 1. Stratified one-dimensional (1-D) seismic and resistivity model structure is shown. For more information see the text
We jointly consider Rayleigh surface waves and electrical resistivities to understand the underground structure underneath the studied area. In fact, there exist many experiments in the literature combining seismic and electrical methods for the underground investigations [22][23][24][25][26]. For the Rayleigh surface waves data we consider combination of single shot and linear geophone spread, i.e., traditional shot-gather. The multichannel analysis of surface waves -MASW [27,28] on such a geophone spread is typically employed to attain the phase velocity dispersion curve. The average seismic velocity structure underneath the geophone spread is obtained by inverting the dispersion curve. The latter velocity structure corresponds to the one-dimensional (1-D) S-wave velocity-depth profile [29]. The MASW method can be applied in two conditions, i.e. passive or active practices in the field. The seismic waves are generated by a weight drop, sledgehammer or vibroseis in the case of active practices. In the case of passive practices (ReMi), the seismic waves generated by cultural or ambient noise (human and industrial activities, construction, machines, cars, waves, winds) recorded by different geophone arrays (i.e., circles, lines, triangles, rectangles) are considered [30,31]. Herein, we consider the phase velocity dispersion curve obtained by the analysis of multiple station (geophone) recordings (i.e., active MASW), and the group velocity dispersion curves acquired by the single-station analysis of waveforms like the approach commonly utilized in the earthquake seismology.
The theoretical modeling that includes both forward and inverse solutions is critical to design an effective field survey with a specific target (e.g., ore body exploration, underground water detection, mapping underground cities, foundation studies) below the surface. The geophysical equipment set, the field layout of the selected equipment and the crew members must be designed well ahead of the deployment to the field for a successful data collection. The design process must include the selection of the right software package that will be employed for the data processing and interpretation. In this respect, we have performed the current theoretical computations where we particularly aim to show our experiences with the data processing steps. By forward solutions we simulate the actual waveforms of Rayleigh surface waves propagating in the Earth by considering a representative model structure. From the computed waveforms are obtained the group and phase velocity dispersion curves. The inversion strategy is based on these dispersion curves to attain the shear-wave velocity variation underneath the studied area. A similar strategy is also followed for the electrical resistivity inversion. On a layered resistivity model representing the Earth structure we compute apparent resistivity pseudo-section, which is then inverted to evaluate the subsurface resistivity distribution.

Rayleigh Surface Waves
The P and S body waves and the Rayleigh and Love surface waves traveling inside the Earth are commonly modelled to study the distribution of compressional and shear velocities within the Earth [32][33][34][35]. The surface waves are primarily sensitive to S-wave (shear) velocities. A sledgehammer (seismic source) vertically hitting the surface mostly creates P waves while the SV waves are created by the wave type conversion from P to SV at the velocity discontinuities (boundaries) in the wave propagation media between the source and receiver. The latter process is accompanied by the constructive interference of P and SV waves multiply refracted and reflected within the Earth layers, which ultimately yield Rayleigh surface waves propagating near the surface. Figure 2 illustrates the shot-gather geometry along with the P-wave source and the vertical geophones. Herein, we do not include the effect of the geophone impulse response. For a particular case, the response curves provided by the manufacturer (i.e., amplitude and phase) can be considered. We can model the actual wave propagation by fitting the observed waveforms to the synthetic waveforms [36]. In this modeling effort, we need the set of model parameters (see Figure 1) and the source information, i.e., the source type (single force, explosion, double-couple), depth, source time history along with the magnitude, and offset (see Figure 2). On the other hand, the surface waves show dispersion. Therefore, we can model the observed dispersion curves (both group and phase velocities) extracted from the waveforms where a match between the observed and the calculated dispersion curves is required. Herein, we employ the second methodology. In the surface wave propagation, there exist several modes (fundamental mode plus higher modes) travelling simultaneously in the media. The fundamental mode -FM is the mode most frequently used in the modeling efforts but to include the higher modes can significantly improve the resolution power of the inversion since the higher modes have greater vertical resolution compared to the FM [37,38].
We obtain the phase velocity dispersion curve from a common-shot gather ( Figure 2). For that we use the technique offered by McMechan and Yedlin (1981) [39]. In this technique, the seismic traces recorded by an array of seismographs are slant stacked in the phase slowness-intercept time domain ( − ) and then transformed into the spectral domain ( − ) by Fourier transform over the intercept time ( ). The fundamental mode plus higher modes superimposed in the recorded waveforms are included in the stacking procedure [40]. The largest values on the phase velocity stack, which is a two-dimensional (2-D) plot as a function of and , are used to select the phase velocity dispersion curve of a specific mode [25,41]. The weighted preconditioned Linear Radon Transform -LRT on the shot gather is also used to acquire the phase velocity dispersion curve. We utilize the LRT alongside the slant stacking -SS. The LRT may occasionally enhance the resolution in the dispersion image by more than 50% (Luo et al., 2008) [42].
Besides the MASW (or SS) method we employ the Multiple Filter Technique -MFT to attain the group velocity dispersion curve [43,44], which represents the average wave propagation from the source to the receiver (geophone). The latter is the single-station approach regularly applied in the earthquake seismology [45][46][47][48][49]. The MFT employs the narrow-band Gaussian filter [− ( − ) 2 2 ⁄ ] on the observed waveforms in the frequency domain. Herein, is the center frequency of the filter. The tradeoff between the frequency and time domains is provided by the filter parameter, i.e., = 35. For a single shot-gather we obtain one phase velocity curve obtained from the SS (or LRT) and multiple group velocity curves obtained from the MFT. The phase velocity curve represents the 1-D average velocity structure beneath the geophone spread. This velocity structure is employed as initial model for the consecutive inversions of the group velocity curves via the iterative least-squares. Therefore, we jointly utilize one phase velocity curve (i.e., MASW) and multiple group velocity curves (i.e., single-station approach) to describe the lateral velocity assembly beneath the geophone spread.

Geoelectrical Method
The geoelectrical properties of the ground depends on several factors, i.e., rock/soil type, intensity and chemical structure of salts, temperature, grain size and distribution of material, pressure and compaction, directional dependence (anisotropy), and moisture level [50]. The electrical conductivity anomalies commencing from variations between neighboring rock/soil units are pursued in the geoelectrical method. Some resistivity values belonging to rocks, minerals, and other Earth materials are listed in Table 1. It is evident that there are great differences between the minimum and maximum resistivity values observed for the rock units or ore minerals. For instance, the quartzite has the resistivity range 500-800000 Ω m. Highly resistive quartzite (i.e., 800000 Ω m) or caliche or laterite layers [51], if located near the surface, would strongly restrict the signal contributions to the deeper layers. On the other hand, the clay has the resistivity range 1-100 Ω m. A near surface clay layer with resistivity (for instance) around 5 Ω m may essentially short-circuit the electrical current permitting very little to penetrate to deeper layers [51]. Actual field conditions may yield in-situ resistivity values outside the ranges presented in Table 1 [53]. The electrical resistivity of sedimentary rocks covers a broad range where the electrical conductivity changes according to the clay and other minerals content, porosity, water saturation. For instance, because of rich clay content shales have low resistivity (5 to 30 Ω m) while limestones, dolomites, and conglomerates show high resistivity (1000 to 100000 Ω me.g., Palacky, 1988) [53]. The geoelectrical method has found several uses in geophysics such as exploring the historic underground cities [12], characterizing the landslide geometry [13], identifying the sedimentary layers [10], determining the bedrock depth [54], studying the sedimentary layers for foundations [14], and searching for groundwater [15].
Various electrode arrangements (or arrays) exist, i.e., Wenner, Schlumberger, Dipole-Dipole, Pole-Pole, Pole-Dipole, Multiple-Gradient [51]. An array is rated effective (i.e., supporting high quality data and well resolved model) if it provides extensive depth penetration with partial electrode expansion, is resilient against the levels of background noise, is less laborious in the field operation and comes with good lateral and vertical resolution. Unfortunately, none of these arrays have these skills all together. For instance, the Wenner array is known effective for the stratified media while the Dipole-Dipole array is used for identifying the lateral changes. The Schlumberger array is frequently favored since only two electrodes are relocated [50]. For a whole length of a traverse, if the field data is collected at several different spacings corresponding to the multiples of a basic electrode spacing, then the effects of lateral resistivity changes can be separated from vertical changes. An array such as Multiple-Gradient supported by substantial computing power and quick field operation may help achieve greatly enhanced images of real resistivity variations. Modern geophysical applications employ Electrical Resistivity Tomography -ERT where the pseudo cross-section measured via a system with multiple electrodes, multicore cables, and multichannel resistivity meter (see Figure 3) is inverted into the subsurface resistivity distribution [55,56]. We employ two programs by Loke (1997) [57] for the forward and inverse modeling of the resistivity data, i.e., RES2DMOD and RES2DINV, respectively. The synthetic model structure, which we consider representing the electrical resistivity-depth profile in the actual Earth, is multilayered over a half-space [58]. The apparent resistivity pseudo-section computed from the RES2DMOD, which is noise free, is then inverted by using the RES2DINV. There are cases that may cause resolution problems in the inversion [12], i.e., coarse electrode spacing, sharp resistivity change from target to host media, and noise. Among the customary data collection arrays (e.g., Wenner, Dipole-Dipole, Pole-Pole, Pole-Dipole, Schlumberger) we utilize the Pole-Pole array to model the current resistivity data. Instead, one may choose to assemble all possible electrode configurations into a single dataset that can be used to enhance the resolve of the resistivity inversion [12][13].

Numerical Results
The actual Earth structure displays three-dimensional heterogeneity. The scale length of these heterogeneities changes from a few tens or hundreds of meters (e.g., inclusions such as sinkhole, cavity, ore body, melt, fluid, coal, salt domes, oil/gas traps in the uppermost crust) observable in the frequency range approximately 1-100 Hz to hundreds of kilometers (i.e., crustal roots, faults, magmatism, mantle upwelling, subduction zones in the crust and upper mantle) in the period range around 5-200 s. Herein, we focus on relatively smooth near surface Earth formations that can be approximated by 1-D (stratified) models. In some regions, the tectonic forces (i.e., compressive, tensional, and shearing forces) severely deform the stratified structure resulting in the thrust-and-fold belts [59] and fault zones where the stratum interfaces are folded and displaced by different faults (i.e., reverse faults, normal faults, or strikeslip faultse.g., see Roche et al., 2020 [60]). Different techniques including three-dimensional seismic, computer and laboratory modelling are needed to image the structural features in geologically complex regions. In our modeling work, we are primarily concerned about how correctly the thickness, depth, and physical parameter (i.e., or ) of a layer are recovered. The electrical resistivity and the group and phase velocities of Rayleigh surface waves are assumed to characterize the Earth structure. Table 2 lists the model parameters defining the lithologic variations in a hard rock environment. The multilayered model made of mainly granite layers (fractured, weathered and water saturated) is assumed to have nine layers with 36-m thickness over the half space (fresh granite). The layer thickness (ℎ), compressional-wave velocity ( ), shearwave velocity ( ), density ( ), shear-wave quality factor ( ) and electrical resistivity ( ) are parameters describing the model structure.

Rayleigh Surface Waves Modeling
We currently utilize the fundamental mode (FM) surface wave propagation to model the underground velocity structure. The FM experiences the interference problem by the higher modes that concurrently travel with the FM. We try to suppress this interference effect in the data processing stage. However, if the FM cannot be differentiated from the higher modes, then the observed FM may be somewhat distinct than the actual curve [63,64]. The seismic surface waves have greater depth penetration with increasing period [65] and the depth penetration is even greater for the higher modes [66]. Therefore, increasing source depth means increasing higher mode influence in the FM dispersion analysis. Earthquakes have source depth from being crustal to upper mantle and in earthquake seismology, mostly shallow earthquakes (i.e., focal depth ≤ 50 km) are employed in the surface wave dispersion analysis so that the interference problem of higher modes could be eased [67].
The software package "Computer Programs in Seismology" by Herrmann (2002) [41] is utilized to compute the dispersion curves and synthetic seismograms, to extract the phase and group velocity dispersion curves, and to invert the dispersion data for the speed-depth profile through the damped least-squares. The slant stacking [39] and multiple filter technique (Dziewonski et al., 1969 [43]; Herrmann, 1973 [44]) are employed to obtain the phase and group velocity curves, respectively. The group and phase velocity dispersion curves corresponding to the multimode Rayleigh surface waves transmitting in the model structure (Table 2) are shown in Figure 4 where the frequency band has the range 2-50 Hz and the fundamental mode -FM plus four higher modes are obvious. The left and right panels show the phase and group velocity curves, respectively. Depending on the underground speed structure there may exist mode kissing effect between different modes (e.g., Gao et al. 2016 [64]), which is not observed in Figure 4 as phase velocity curves are well separated. In terms of group velocities, there exist some higher modes mixing at high frequencies (>30 Hz) and at low group velocities (<900 m/ssee the ellipse in the right panel). Note that the FM is free of mode mixing from the higher modes. However, this may not always be true since some higher mode energy (particularly first higher mode) may converge to the FM in which case the interpreter should be extra cautious against mode confusion while selecting the FM curve. Figure 5 displays the synthetic seismograms travelling in the considered model structure ( Table 2). In the computation, the wavefield is sampled at 0.002 s and is recorded by 48 vertical geophones with 2-m constant spacing. The minimum and maximum offsets are 2-m and 96-m, respectively.  Table 2 are shown. The fundamental mode -FM is labeled n=0 and then higher modes are labeled with increasing numbers, i.e., n=1, n=2, n=3, and so on. The phase and group velocities are shown in the range 300-1700 m/s. Table 2. See text for more explanation

Figure 5. Seismic traces showing P-SV body waves and Rayleigh surface waves corresponding to the model structure in
We have applied the Slant Stacking -SS to the synthetic seismograms shown in Figure 5. Figure 6 (left panel) shows the resulting fundamental mode -FM phase velocity curve specified by a white line connecting black dots that determine the maximum stack (red color area) on the dispersion image. The right panel in Figure 6 displays the consequence of Linear Radon Transform -LRT with weighting and preconditioning. Superimposed on the dispersion image (blue lines, right panel) are shown the theoretical multi-mode phase velocity curves corresponding to the model structure in Table 2. For frequencies less than 8 Hz the FM phase velocities estimated by the SS (left panel) or LRT (right panel) show some uncertainty, otherwise the solid circles (LRT) and the blue line (theoretical) shows good fit. The geophone spread in Figure 5 corresponds to a lateral extent of 94-m where the first geophone is located at 2-m, and the last geophone at 96-m. The FM phase velocity curve displayed in Figure 6 characterizes the (average if heterogeneous) velocity structure beneath the geophone arrangement.
We invert the FM phase velocity curve in Figure 6 to retrieve the 1-D shear-wave speed structure (Table 2) Figure 7 is performed under noise free conditions. Note that because of observational uncertainties at periods longer than 0.1 s the dispersion fit between the inverted (red color) and theoretical (open circles) is relatively weak. The inverted (red color) model structure is employed as initial in the subsequent 1-D inversions for which we employ the group velocity curves. Table 2

The horizontal axis shows the period (s)
We find the fundamental mode -FM phase velocity curve ( Figure 6) by applying the multichannel analysis of surface waves (MASW) for which we use two equivalent methods (i.e., Slant Stacking -SS and Linear Radon Transform -LRT). The next step is to apply the traditional single-station method to obtain the FM group velocity curves related to the pathways between the source and the multiple receivers (geophones). This application needs source information, i.e., source location and timing where a coordinate system in the field is used to calculate the source and geophone locations (see Figure 2) and the trigger cable precisely transmits the source timing (or zero time) to the seismograph. The constructive interference mechanism mentioned above determines the frequency band of propagating surface waves. The short wavelength (or high frequency) surface waves are first created in the near field and then the frequency band broadens to long wavelengths (or low frequencies) in the far field, i.e., progressively away from the source. The depth penetration is determined by the wavelength, i.e., surface waves reach deeper with growing wavelength.
We select a distance range to utilize the single-station surface wave recordings from 30-m to 90-m, i.e., these geophones numbered from 15 to 45 (see the stars in Figure 5). This distance range is effective to determine the speed structure in the 36-m thick model (Table 2) above the half-space. A range of wavelengths is practically utilized to invert the underground speed structure. At distances shorter than 30-m these long wavelength surface waves are not constructively created obstructing the determination of the deep structure. At distances longer than 90-m these short wavelength surface waves are inelastically attenuated and then the near surface speed structure is not resolved well. For target depths deeper than 36-m the interpreter may need to consider these receivers (geophones) at distances longer than 90-m to include the longer wavelengths in the dispersion analysis. Figure 8 displays the outcome of multiple filter technique -MFT applied to the geophone recording at 60-m distance (i.e., receiver #30 shown by a tick mark in Figure 5). The technical details regarding the MFT are provided in Herrmann (2002) [41]. In Figure 8, the white line linking the solid squares indicating the maximum amplitudes (i.e., middle of the red color area) on the dispersion image illustrates the FM group velocity curve. To enhance the relatively weak amplitudes in the waveform, we treat the surface wave spectral amplitudes | ( )| in the MFT as follows.  In Figure 9, we invert the FM group velocity curve (Figure 8) for the 1-D shear-wave velocity structure corresponding to the wave propagation between the source and geophone #30. The initial model accepted from Figure  7 is represented by the blue dashed line and the inversion result by the red line. In Figure 9, the inverted (red line) and the initial (blue dashed line) speed structures closely fit. This is expected since the theoretical model structure in Table  2 is one-dimensional (1-D). However, the two speed structures are still somewhat separate from each other particularly below the 32-m depth. This is partly due to the group and phase speeds having distinct resolution control varying with depth [70]. There are some mismatches between the theoretical (observed) group velocities shown by the triangles and the inverted group velocity curve indicated by the red line. The uncertainties at particularly low frequencies (or long periodssee Figures 6 and 8) are due to the MFT resulting the observed mismatches in Figure 9.

The horizontal axis shows the period (s)
The single-station inversion strategy is established as in Figure 9. Likewise, we move to inverting the other geophone recordings in the array, i.e., from geophone #15 to geophone #45, which yields 31 FM group velocity curves like the traditional single-station method. The procedure outlined in Equations 2 to 5 is used to transfer these singlestation group velocity curves into grid-based dispersion curves with constant 2-m grid spacing, which is equal to the geophone spacing. Currently the linear system in Equations 2 to 5 has 31 columns consistent with the number of grid points and 31 rows corresponding to the number of single-station (observed) group velocity curves. The respective data processing system developed in Equations 2 to 5 is illustrated in Figure 10. The number of grid points and geophones may vary with respect to the specific problem under consideration.
We continually solve the system in Equations 2 to 5 for each surface wave period so that the period dependent dispersion curves at each grid point are constructed.
where m is the number of geophones and grid points. with dimension ( ) and lower triangular form is the system matrix where the elements ( ) of matrix is one for ≥ and zero for < . ̅ represents the single-station (observed) group slowness, and stands for the group slowness at a grid point attained by solving the linear system above, which can be accomplished by the damped least squares method, i.e., ( + ) =̃ where is the damping parameter used to restrain the noise effect, is the unknown slowness vector, and ̃= 1 + 2 is the known slowness vector corresponding to the righthand side in Equation 2. The average (group) slowness related to the propagation media between the source and the first geophone (i.e., currently geophone #15) is characterized by 0 , which is calculated from the 1-D velocity-depth profile acquired from the inversion of the phase velocity curve (see Figure 7). The 1-D speed structure acquired from the inversion of the group velocity curve at geophone #15 can be alternatively used to calculate 0 . Note that to decrease the geophone spacing (i.e., ∆ ) would increase the horizontal resolve in Figure 10. In Figure 11, using a fine color scale we show the cross-section representing the shear-wave speeds (m/s) of the stratified 1-D model structure ( Table 2). The layer interfaces along with the speeds in each layer are overlaid on the cross section, which is achieved by the application of Equations (2)(3)(4)(5) to the geophone recordings (see the stars in Figure 5). The model structure in Table 2 includes low and high velocity zones (i.e., LVZ and HVZsee Figure 11), which are the features satisfactorily recovered through the inversion procedure defined in Equations (2)(3)(4)(5). However, there are still some resolution losses regarding the interface depth and sharpness. For instance, the near surface layer interface at 4-m depth is vaguely inverted although this depth range is probed by the high frequency (or short wavelength) FM surface wave energy. The high velocity zone (i.e., HVZ) is correctly identified by two layers in the depth range 8-16 m, but the deeper low velocity zone (i.e., LVZ) lacks some precision since the shear-wave speed in the seventh layer underlying the LVZ is underestimated, i.e., ~975 m/s instead of the actual 1000 m/s. The last interface topping the half-space is placed 2-m shallower than the actual depth at 36-m. Because of dispersion ambiguities at low frequencies (see Figures 6 and 8) the shear-wave speed in the half-space is overestimated at around 2150 m/s against the actual speed at 1600 m/s (see Table 2).  Table 2 is shown

Resistivity Modeling
We employ the finite differences source code package by Loke (2014) [71] to calculate the theoretical apparent resistivity pseudo-section (forward modeling) of the stratified resistivity structure in Table 2. In the finite difference method, we utilize grid cells equal to one half the electrode spacing. The profile expansion is 129 m achieved by 44 electrodes arranged at 3-m electrode spacing along with Pole-Pole array. Noise free conditions are considered in the numerical calculations (apparent resistivity values).
There exist three practices (i.e., electrical horizontal profiling, electrical vertical profiling, and electrical vertical sounding) applied in the field to explore the subsurface resistivity structure [50]. The near surface properties of the underground are explored using the horizontal profiling. The latter approach is applied using a fixed electrode spacing and the electrodes are shifted in the field to complete the apparent resistivity survey. The horizontal profiling provides only restricted stratigraphic image in the vertical direction, which can be overcome by the vertical profiling where an arrangement of inner electrode and electrode spacings are selected. In the vertical profiling, the electrode spacing is repeatedly increased and at each step deeper sections of the subsurface are measured to build the corresponding apparent resistivity pseudo-section. In case of the uniform flat layers, one may prefer the vertical electrical sounding (VES) where the electrode deployments in the field include regular increments of depth penetration.
We utilize the source code package by Loke (2002Loke ( , 2004 [72,73] to calculate the actual resistivity structure by inverting the apparent resistivity pseudo-section shown in Figure 12. The inversion algorithm uses the least-squares method along with smoothing in two-dimensional media. The number of numerical computations in the modeling is reduced by using a quasi-Newtonian procedure [74]. The layer thicknesses and resistivities particularly near the surface are effectively inverted. The layer interfaces along with the resistivities in each layer ( Table 2) are overlaid on the cross section in Figure 13. The change of the inverted resistivities with depth is displayed using a rainbow color scale.  Table 2 is shown in terms of the true (inverted) resistivity values Figure 13 shows the inversion result (cross section) for the resistivity model in Table 2. Resistivities (i.e., 150-250 Ω m) in the first and second layers are resolved on the cross section, but these layers (i.e., seventh, eighth, and ninth layers) lack some precision since the inverted resistivities are lower compared to the actual resistivity values in these layers (i.e., 1400-2000 Ω m). The inversion results show that high resistivities (i.e., 1200-1400 Ω m) in the third and fourth layers and low resistivity (i.e., 300 Ω m) in the fifth and sixth layers are not visible on the resistivity cross section. This happens because the low resistivity layers (i.e., fifth and sixth layers) are sandwiched between the high resistivity layers (i.e., fourth and seventh layers) in which case the electric currents tend to accumulate in the middle low resistivity layer. In the opposite case, the high resistivity layers (i.e., third and fourth layers) are sandwiched between the low resistivity layers (i.e., second and fifth layers). The electrical currents are now repelled by the high resistivity zone. The latter two cases result inhomogeneous electrical current distribution in the subsurface and subsequently the inverted resistivity model poorly resolves some depth sections.
The facts about Figure 13 emphasize the importance of joint interpretation of the underground structure. For instance, the magnetotelluric -MT method also shows sensitivity to the electrical properties of the subsurface. However, this sensitivity is mostly related to the deeper parts of the Earth while only limited MT sensitivity exists for the near surface. Therefore, these two methods (i.e., electrical resistivity good for near surface and magnetotelluric good for deeper parts) can be joined together for broader vertical resolution of the underground electrical structure.

Discussions and Conclusion
We have chosen a one-dimensional model structure ( Table 2) to perform some numerical experiments regarding seismic surface waves and electrical resistivity. The surface wave part includes Rayleigh surface waves, which consists of fundamental mode -FM plus higher modes. We have currently considered only FM data while the higher mode energy is restricted from the waveforms by placing the source on the free surface. The higher modes alongside the FM turn into detectable on the phase velocity dispersion image if the seismic source is placed deeper below the surface, e.g., 5-m. Figure 14 displays two examples of the phase velocity dispersion image corresponding to source depths at 5-m (left panel) and 10-m (right panel). The higher modes are visible at high frequencies greater than ~50 Hz. In case of the source depth at 5-m, only two higher modes (first and second higher) come into the display (see left panel). The source depth at 10-m produces extra two higher modes (third and fourth highersee right panel). In the meantime, the frequency band, in which the first and second higher modes are visible, shifts to lower frequencies. The latter is also valid for the FM. The higher modes have greater depth penetration and offer better horizontal and vertical resolution of the subsurface because of their short wavelengths.
We primarily focus on the near surface Earth structures with multiple flat-lying layers. The present approach can be extended to include multi-dimensional heterogeneities (or inclusions such as melt, ore body, sink hole, cavity, salt, fluid) embedded in the stratified Earth. The resolution power of the fundamental mode -FM can be increased by additionally considering the higher modes in the modeling efforts. However, to excite the higher mode Rayleigh surface waves is expensive since the P-wave source (see Figure 2) must be placed at some depth below the free surface level. Digging a trench or drilling a hole as deep as 5-m could be a solution. In the current study, we employ the fundamental mode dispersion curve (both group and phase velocity). The respective dispersion curve needs to be isolated on the dispersion image. If different modes interfere with each other, then the mode isolation becomes difficult [75]. For instance, higher mode group velocities in Figure 4 interfere with each other as indicated by an ellipse. The mode interference event is relatively modest for the phase velocity dispersion curves where mode kissing event mostly occurs [64]. The slant stacking or linear radon transform successfully isolates the FM phase velocity curve (e.g., see Figures 6 and 14), which could be utilized to assess where the multi-mode group velocity curves are expected on the relevant dispersion image. On the other hand, the full waveform inversion -FWI, which does not involve mode isolation, can be alternatively used, but has some drawbacks, i.e., requires source information and excessive computational time (CPU). When the wave transmission media comprises of powerful heterogeneities, the FWI may be unsatisfactory [63].
There exist electrode arrangements (e.g., Dipole-Dipole, Pole-Pole, Schlumberger, Wenner) traditionally used in the resistivity field surveys. Recent resistivity surveys employ more advanced acquisition systems, which are equipped with multiple electrodes, multicore cables, and multichannel resistivity meter [12,13]. The new system can collect thousands of data in a manner of tens of minutes, which brings great sensitivity to the resolution of the subsurface electrical structure. All electrodes except the current electrodes in the multiple electrode array systematically measure the relevant voltages once the current is injected into the ground. One can still apply one of the traditional electrode arrays to the collected data by selecting the appropriate data set for the intended electrode array. Herein, the setting with 44 electrodes, 3-m electrode spacing, and 129-m expansion is used. Both lateral and vertical electrical resolution can be increased by increasing the number of electrodes to (for instance) 96, decreasing the electrode spacing to 2-m, and increasing the expansion to hundreds of meters. We have limited computational power and our resistivity modeling is restricted to what is presently demonstrated.
Modern geophysical applications develop in a way that different geophysical methods (e.g., resistivity, gravity, seismic) with different sensitivity (i.e., conductivity, density, velocity, respectively) are applied in an integrated manner in the field. Therefore, the subsurface lithologies are better delineated. For instance, the occurrence of fluids in the subsurface may be indicated by high electrical conductivity and low seismic velocity [76]. On the other hand, the existence of ore deposits may be revealed by low electrical resistivity and high seismic velocity while the presence of voids may be indicated by high electrical resistivity and low seismic velocity. Therefore, multiple geophysical methods rather than one of these methods may be more beneficial if applied jointly.

Data Availability Statement
The data presented in this study are available in article.

Funding
The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of Competing Interest
The authors declare that there is no conflict of interests regarding the publication of this manuscript. In addition, the ethical issues, including plagiarism, informed consent, misconduct, data fabrication and/or falsification, double publication and/or submission, and redundancies have been completely observed by the authors.