Interpretation Methods for Seismic Downhole Test in Inclined Boreholes

Geotechnical investigations often involve inclined boreholes, which can be used for downhole (DH) seismic surveys. However, as there is no interpretation method for downhole tests in inclined boreholes (IDH), this study proposes alternative interpretation methods based on the direct method (DM), interval method (IM), modified interval method (MIM), and refracted ray path method (RRM). We have named the proposed methods, adding an I to the original name to indicate that they are performed in an inclined well, i.e., DMI, IMI, MIMI, and RRMI. To analyze the applicability of the proposed methods, eight simple models with horizontal layers and four 2D models were used to obtain the P-and S-wave velocity profiles. Among all the proposed methods, the RRMI method showed the best fit between the calculated S-wave velocity (Vs) profile and the real models, providing good reliability. To test the equations and hypotheses, new interpretation steps were developed based on Snell's law and a modification of the numerical bisection method, which showed that the error increased slightly as the dip angle of the well decreased. The next step was to test the accuracy of the RRMI method in the field and develop downhole test processing software for vertical and inclined boreholes.


Introduction
Seismic geophysical tests are widely used to determine P-and S-wave propagation velocity profiles and soil properties that are essential for evaluating the dynamic behavior of soil deposits, the geotechnical design of foundations, and the seismic response analysis of geotechnical structures [1][2][3].These tests require the installation of seismic sensors (geophones or receivers) and an energy source for generating seismic waves (a hammer or explosive charge).The receiver can be placed on the ground surface or borehole wall.The most common seismic geophysical tests performed at the ground surface are seismic refraction and surface wave measurements in multichannel arrays (MASW, SPAC), whereas the tests performed in boreholes are the downhole (DH), cross hole, and uphole tests; the DH test is one of the most commonly used.Conventionally, DH tests are performed in a vertical borehole, where the P-and S-wave travel times are measured from a surface seismic source to a geophone located at different depths within the borehole [4].Alternatively, DH tests can be performed using a seismic cone penetrometer (SCPT) rig.In this case, the geophone was placed at the tip of the penetrometer, and the test was performed as the penetrometer vertically reached different depths [5,6].Figure 1 shows a schematic diagram of the execution of the DH test in a vertical borehole.
Currently, seismic geophysical tests are frequently performed as part of geotechnical investigations and are an essential complement to direct exploration techniques [7][8][9][10], particularly in countries with high seismic activity.In these environments, it is necessary to determine the dynamic behavior of the soil to evaluate the problems of seismic amplification and liquefaction of soils or to analyze the seismic response of structures.Detailed geotechnical exploration programs developed for these purposes, mainly in the mining industry, frequently include inclined boreholes, many of which avoid intangible zones or risk zones, as well as determining the mechanical properties of a specific zone (rock mass quality, permeability, etc.) [11][12][13].However, it is impossible to perform DH tests in these boreholes because there are no interpretation methods for this type of borehole [14].

Figure 1. Field scheme of the conventional Downhole
This research proposes methods for interpreting downhole tests in inclined boreholes (IDH) by reformulating existing interpretation methods and using synthetic data on elastic wave propagation in a stratified soil model (Figure 2).Wave propagation data were obtained for a simple horizontal layered soil model using Snell's law and a complex layered soil model using the finite difference method (FDM).

Research Methodology
To develop interpretation methods for IDH tests, methods used in practice have been analyzed and modified owing to the simplicity of their algorithms [15,16], including the Direct Method (DM) [17], Interval Method (IM) [18], Modified Interval Method (MIM) [18], and refracted ray-path method (RRM) [19].The methodology used to achieve the objectives of this study is described in the flowchart shown in Figure 3.

Define conditions and limitations
Generate sinthetic data based on Snell s law and elastic wave propagation modelling by FD Evaluate the formulated methods for inclined borehole using sinthetic data

Analysis and discussions of results
Analyze the performance and limitations of the proposed methods

Review of Developed Methods
The DM, IM, and MIM methods assume that the seismic ray path is in a straight line through the bedded medium and that the horizontal distance from the seismic source to the receivers (along the borehole) is much smaller than the depth of each receiver.However, seismic waves generated at the surface travel through soil strata with different stiffnesses and therefore refract at each interface according to Snell's law and Huygens' principle [20].Therefore, many researchers recommend caution when using these methods in heterogeneous media or with a higher impedance contrast [21].The RRM method, unlike the previous methods, is based on Snell's law and considers that seismic rays change direction owing to the change in incidence angles in each layer [17]; Consequently, it considers the travel times of each layer (separation between receivers) along the ray path to the receivers, which greatly reduces the error in the results of P-and S-wave velocities.Because these methods are commonly used in practice [18,19], the evaluation and modification of the interpretation of IDH test data were performed in this research.

Initial Formulation and Conditions
The following conditions and limitations of the physical properties of the materials and laws of seismic ray propagation in stratified media were considered for the formulation of the interpretation methods of the IDH test:  The horizontal distance from the seismic source at the ground surface to the receivers (along the borehole) was less than the depth of each receiver.
 The investigated subsoil consisted of homogeneous horizontal strata without lateral variations in elastic properties.
 Each soil layer has isotropic elastic properties; thus, in each homogeneous layer, the seismic wavefronts are circular, and the trajectories of the seismic rays are always correctly distributed.
With these considerations, the DM, IM, MIM, and RRM interpretation methods were reformulated for processing the IDH test data.The new names for these methods are DMI, IMI, MIMI, and RRI, where the added I indicates that the test was performed in an inclined well.

Direct Method for Inclined Borehole (DMI)
The travel time of the seismic wave from the source to the receiver along the inclined borehole was obtained like that in a conventional DH test, assuming that the ray path is rectilinear.This travel time was corrected by assuming that the horizontal distance from the source to the receiver was less than the depth of each receiver, using Equations 1 and 2, as shown in Figure 4 (left).
where tc is the corrected travel time, D is the depth of the receiver along the inclined borehole, t is the travel time of the wave from the source to the receiver, R is the distance from the source to the receiver,  is the angle of inclination of the borehole, S is the distance from the source to the wellhead, and Z is the vertical depth from the ground surface to the receiver.

Interval Method for Inclined Borehole (IMI)
The travel time interval between each pair of receivers was obtained in the same manner as in conventional DH tests by subtracting the travel time of the upper receiver from that of the lower receiver, as shown in Figure 5. Considering that there is no lateral variation in the material properties and that this method is based only on the travel time and total distance from the source to each receiver, the wave velocity in the interval between the upper and lower receivers can be calculated using Equation 4: The straight-line ray path R was calculated in the same manner as in the direct method using Equation 2, where V corresponds to the wave velocity in the interval between the upper and lower receivers.

Modified Interval Method for Inclined Borehole (MIMI)
Similar to previous methods, this method is based on the rectilinear path of seismic rays and assumes that the horizontal layers have the same thickness, proportional to the interval between the receivers.The propagation time of the seismic wave from the source to each receiver depends on its velocity in each soil layer, which is not necessarily constant, as shown in Figure 6.Consequently, the time it takes for the waves to reach the upper and lower receivers is equal to the sum of all the time intervals that the wave takes to pass through each soil layer and is determined by Equations 5 and 6.The distances traveled by the seismic waves in each layer were calculated using Equation 7. Finally, the wave propagation velocity for each interval was calculated using Equation 8.
where Vi is the seismic wave velocity for the i th layer, Zi,l and Zi,u the vertical depths in the i th layer of the lower and upper receivers, respectively, at this position, Ti,u and Ti,l are the total travel times in the i th layer from the source to the upper and lower receivers, respectively, Lij,l and Lij,u: Lengths of the ray path segment in the j th layer traversed by the seismic ray directed towards the i th layer at the lower and upper receivers, respectively, Ri,l is Total length of the ray path from the source to the lower receiver for a seismic ray directed towards the i th layer.This value can be calculated using Equation 2. This method is similar to the MIMI method; the only difference is that the seismic ray path considers wave refraction according to Snell's law.Consequently, the ray trajectory must satisfy the Equations 9 to 13: where ei is the thickness of each layer, and θij is the angle of incidence of the seismic ray for each layer.Equation 9 is based on Snell's law, which considers the energy dispersion when a plane wave is incident on the contrast of the impedances.Equation 10represents the geometric equality between the horizontal distance traveled by the seismic ray to the ith receiver and the horizontal distance from the source to the receiver, as shown in Figure 7.In the first layer, the interval velocity (V1) is calculated by assuming a rectilinear ray path with the first receiver.The velocity calculation of the subsequent layers requires an initial velocity model, which can be determined using the DMI, IMI, and MIMI methods or from a random velocity model.Then, an iterative process is performed, as shown in the processing scheme of the RRMI method in Figure 8 and detailed in the following steps: 1.The value of θi1 is randomly defined.9was used to calculate the value of θi2 (V2 was obtained from the initial velocity model):

Equation
3. The error between both terms of Equation 10 is obtained.11was used to calculate the lengths of the seismic rays in each layer (Li).

Equation
5. Using Equations 12 and 13, the new value of V2 was calculated, and the velocity model was updated.
6. Using V1 and the new V2, a new value of θi1 is calculated, and steps 2-6 are repeated until the error in step 3 approaches zero.
7. The next record is then processed by repeating the same process.
To automate this process and find the correct θi1, a modification has been made to the "bisection method," which was introduced shortly after the intermediate value theorem proposed by Bernardo Bolzano in 1817 [22].The bisection method is a simple but robust root-finding algorithm that repeatedly divides an interval into two equal parts and then selects the subinterval in which the root exists [23].According to the bisection method theory, the objective function f(x) must be continuous and will be equal to zero only if it has a solution [23], while the objective function θi1 is continuous in the range of [0° to 90°] According to Fermat and Huygens-Fresnel principles, there exists a single seismic ray path with a specific θi1 that arrives first at the receiver in an elastic-homogeneous medium [24].The iterative process proposed for the RRMI is shown in Figure 8.

Figure 8. Processing scheme of Refracted Ray-path Method in Inclined Borehole (RRMI)
The bisection method allows high accuracy in the iterative process because the only existing solution is within the search range of this method.As the iterations progressed, the method approached the solution with fast convergence, as shown in Figure 9.

Comparison of the Proposed Methods Using Synthetic Models
To verify the reliability of the proposed methods, synthetic velocity models were used based on the path of seismic rays from the source to the location of each receiver along an inclined borehole according to Snell's law.Two-and threelayer soil models were considered for wells with dip angles of 30°, 45°, and 60°.The separation between receivers was set to be 1.0 m along the borehole, and the horizontal distance from the source to the wellhead was 3.0 m, as shown in Figure 10.These models extracted the travel time of seismic rays from the source to each receiver using Snell's law.Subsequently, the shear-wave velocity profile along the borehole was obtained using the DMI, IMI, MIMI, and RRMI methods.For each of these velocity models, the travel time of the seismic waves to each receiver was calculated using Snell's law, as listed in Table 1.Subsequently, the proposed DMI, IMI, MIMI, and RRMI methods were applied to determine the velocity profiles of the shear waves (Vs) using the arrival times of the waves obtained for each model, as shown in Figures 11 to 13.The results show that none of the methods have uncertainty when calculating the velocity of the waves in the first layer because this is homogeneous and the ray paths are rectilinear.Different degrees of uncertainty were found in calculating the velocity of the waves in the second layer, depending on the method used.For the case in which the velocity of the second layer was greater than that of the first, the maximum average errors found in the calculation of the velocity of the second layer were 27.0% with the DMI method, 19.2% with the IMI method, 11.0% with the MIMI method, with an exceptional average error of 113% in one case, and 5.7% for the RRMI method, in which 83% of the cases analyzed with this last method presented average errors lower than 0.1%.In the case where the wave velocity of the second layer is lower than that of the first layer (Vupper > Vlower), the velocities of the second layer calculated by all the evaluated methods are much more accurate, possibly because of the precision of the determination of the first incident ray during the iterative process.It is also observed that as the borehole inclination increases, the borehole becomes more vertical, and the errors in all the methods used decrease.

Three-Layer Model
In the case of the two-dimensional soil model consisting of three layers, a well with an inclined depth of 15 m was proposed, and four models of shear wave propagation velocities were defined: Vs1 = 300 m/s, Vs2 = 500 m/s, and Vs3 = 800 m/s for the first model; Vs1 = 300 m/s, Vs2 = 800 m/s, and Vs3 = 1200 m/s for the second model; Vs1 = 800 m/s, Vs2 = 500 m/s, and Vs3 = 300 m/s for the third model; and Vs1 = 800 m/s, Vs2 = 300 m/s, and Vs3 = 500 m/s for the fourth model.In each model, the travel time of the seismic waves to each receiver was calculated using Snell's law, as shown in Table 2. Subsequently, the proposed DMI, IMI, MIMI, and RRMI methods were applied to determine the velocity profiles of the shear waves (Vs) using the arrival times of the waves obtained for each model, as shown in Figures 14 to 16.The results show that for soil profiles where the velocity increases with depth, good accuracy is obtained between the synthetic velocity model and that calculated by the RRMI method, and to a lesser extent with those obtained by the IMI and DMI methods.The MIMI method did not show a good fit; it became unstable when there was a strong contrast between the velocities of adjacent layers.However, all methods showed a better fit to the synthetic model in soil profiles with decreasing velocities with depth.The maximum average errors were 25.2%, 14.0%, 35.6%, and 0.8% for the DMI, IMI, MIMI, and RRMI methods.It is noteworthy that for 83% of the models analyzed, the maximum average error found with the RRMI method is 0.1%, which is the method that gives the best results in the case of inclined wells.As in the case of the two-dimensional soil model with two layers, all methods show lower average errors as the angle of inclination approaches the vertical direction, as the horizontal distance from the source to each receiver decreases, and the ray path, therefore, has a greater linearity, approaching that assumed by the interpretation methods.

Discussion
The results obtained for the two-dimensional soil models with two and three layers show that the average error of the velocity models obtained using the proposed processing methods is strongly influenced by the dip angle of the borehole, as shown in Figure 17 (note that the y-axis is on a logarithmic scale).This is because the smaller the dip angle of the borehole, the more refraction the seismic waves suffer between layers; thus, the DMI, IMI, and MIMI methods, which assume straight ray paths, have higher errors than the RRMI method, which considers the refraction of the seismic waves in each soil layer.This condition is more pronounced when the propagation velocity of the seismic waves increases with depth, whereas when the velocity decreases with depth, all the proposed methods converge rapidly.The errors decrease as the borehole becomes more vertical, as shown in Figure 18.The best approximation of the seismic profiles was obtained with the proposed RRMI method, which yielded negligible errors; therefore, this method is the most recommended for interpreting the seismic profiles of an IDH test.The iterative algorithm used in this method to search for the angle of incidence is simple.It requires a shorter search time because the optimal angle is generally found near the middle (between 0° and 90°) and not at one of the extremes (near 90°).However, it was found that when applying the RRMI method, the choice of the initial velocity model was a critical issue for improving the results and reducing uncertainty, as shown in Figure 19, where the initial velocity models obtained by the DMI, IMI, MIMI, and a default model were used.

Figure 19. Vs models calculated by RRMI method using the DMI, IMI, MIMI and a default model as the initial model
However, it has been found that in the models analyzed, especially those in which the velocity increases with depth, the MIMI method presents significant errors and even becomes unstable when there is a strong contrast in velocity between soil layers; therefore, its use is not recommended for processing the IDH test.

Lateral Variation of Shear Wave Velocity and Wave Modeling
Another critical point for evaluating and verifying the applicability of the proposed interpretation methods for IDH tests is the effect of lateral variations in soil properties.Synthetic records of the propagation of elastic seismic waves in a heterogeneous medium were generated using the finite-difference method (FDM).The first-arrivals were extracted from these records, and the travel-time of the waves from the source to the receivers was obtained.
FDM has become one of the main tools for solving seismic wave propagation problems [25].This method is based on the representation of a defined continuous function τ(x, y, z, t) and its derivatives in terms of an approximation τ(tau) defined at certain discrete points called grid points [26].
The plane-wave Equation 14 is based on an eikonal wave, a first-order nonlinear partial differential equation [27].To solve the wave propagation in 2D models and discretize the number of time samples and sampling intervals ( and , respectively) during the modeling by FDM, stability criteria and first-order equations according to Kramer [19] were considered.In addition, the Ricker wavelet (Equation 15), which is often used in seismic analyses [28,29], was used to generate the first pulse of the source, where f is the peak frequency.
Three seismic profiles of the 2D-MASW tests provided by ZER GEOSYSTEM PERU SAC were selected to obtain representative cases of two-dimensional shear wave profiles.The 2D-MASW is a noninvasive geophysical test based on surface wave propagation, which allows us to determine the lateral variation in shear wave velocity [30,31].Three two-dimensional section models (M-01, M-02, and M-03) were extracted from these profiles, each with a different lateral variation in the shear wave velocities (magenta boxes in Figure 20-a).The size of these models was 40 × 40 m, gridded with 0.1-meter elements.The wave generation source was located 3 m from the top of a hypothetical well with a dip angle of 37°, and receivers were spaced 1 m along the well, as shown in Figure 20-b From the synthetic wave propagation modeling, the arrival time of the wave at each receiver was recorded for each analyzed model (Figure 21, upper part).Subsequently, using the proposed method, the velocities of the shear waves along the inclined borehole were calculated to obtain a 1D velocity profile, as shown in the middle and bottom of Figure 21.
From the results, it was observed that the Vs profiles obtained by the RRMI method fit better to the theoretical models extracted from the 2D MASW tests, with errors of 5.1% and 5.0% for the M-01 and M-02 models, respectively.In contrast, the lower part of the Vs profile of the M-03 model did not present a good fit, and the error increased to 7.4%.This increase in error may be due to the strong lateral variation in the shear wave velocities presented by this twodimensional model.The IMI method also presented relatively low errors, with values of 5.3%, 11.8%, and 8.7% for the M-01, M-02, and M-03 models, respectively.By contrast, the DMI and MIMI methods present higher errors in estimating the Vs profile; therefore, their application in these cases is not recommended.As a second case study to test the proposed methods, a 2D synthetic model consisting of complex soil layers (recent soil layers, dipping layers, fault zones, and angular unconformities) was considered (Figure 22-a  The first arrival time was manually selected from the wave records at each receiver, and the S-wave velocities were calculated using the proposed method, as shown in Figure 22-d.The DMI method could not be used in this case because the time-distance curve had a sinusoidal shape, which made it difficult to identify the layers.The IMI, MIMI, and RRMI methods present errors of 30.9%, 12.8%, and 18.7%, respectively, which are higher than those found in the three models analyzed in the first case, probably because of the stratigraphic complexity of the model.
From the results obtained for the analyzed cases, it was observed that the shear wave velocity profiles obtained by the proposed RRMI method fit better with the synthetic models, both for simple and complex stratigraphic models.The low uncertainty presented by this method is due to several factors, such as the principle used, the sophisticated algorithm, the numerical method of the iterative process, and the initial model used.This last factor is one of the most relevant factors in processing the RRMI method; therefore, it is recommended to use an initial velocity model obtained by the DMI, IMI, or MIMI.

Conclusions, Recommendations, and Future Steps
 Methods for interpreting inclined downhole seismic surveys (IDH) were proposed by modifying the existing DM, IM, MIM, and RRM methods for vertical boreholes.
 Among the proposed interpretation methods for inclined downhole tests (IDH), the DMI, IMI, and RRMI methods generally show good reliability in obtaining shear wave velocity profiles in soils with horizontal layers, with the RRMI method having the lowest error.The MIMI method has some inconsistencies when the wave velocity in the soil layers increases with depth and becomes unstable when the wave velocity contrast between two consecutive layers is large.
 The borehole inclination affected the accuracy of the results for all methods; the more vertical the borehole, the more accurate the calculations.This is because the verticality of the borehole reduces the difference between the angles of incidence and refraction in each layer; thus, the assumption made by the methods in which the wave travels in a straight line is more accurate in this case.
 The FDM-2D modeling results for the different cases show that the lateral variation of the soil layers (lateral changes in thickness and velocity) and the spatial configuration of the refractory interfaces (dip of layers, erosion surfaces, fault zones, angular unconformities, etc.) strongly affect the calculation of wave velocities using the proposed methods.
 According to the analyses and results obtained, the RRMI method is the best option for interpreting the seismic data of the inclined downhole test (IDH) because, in these wells, the horizontal distance from the source to each receiver causes the seismic wave to have greater refraction at the interface of each layer.This effect became more significant as the receiver depth increased.
 The selection of the initial velocity model for the RRMI method significantly influences the results; therefore, good criteria are required for the selection of the method to be used to generate the initial model, either DMI, IMI, or MIMI.Therefore, for optimal processing of shear wave velocities, it is necessary to calibrate the model with a borehole log (material contact, material type, etc.).Using a random-velocity model as the initial model for the RRMI method is not recommended.
 As this research was conducted using synthetic data, the next step was to perform tests in inclined boreholes to confirm the accuracy of the proposed RRMI method, which will obtain new parameters that may not have been considered in this research.
 The results of the RRMI method were favorable for the vast majority of the evaluated synthetic models (horizontally layered and laterally varying stratified soil profiles).Therefore, based on the results of this research, we plan to develop seismic processing software for IDH tests.

Data Availability Statement
The data presented in this study are available on request from the corresponding author, as they belong to the consulting services of the ZER Geosystem Peru S.A.C.

Funding
This work was supported by ZER Geosystem Peru S.A.C.

Figure 2 .
Figure 2. Scheme for a downhole seismic survey of an inclined borehole

Figure 4 .
Figure 4. Direct Method in an inclined well.Left: Scheme of the test execution.Right: Corrected time plot to estimate the wave velocity from the slope of the lines Next, tc versus D was plotted for each receiver along the borehole, and the S-wave velocity of the layers was obtained from the slope of the line generated from the data fit using Equation 3, as shown in Figure 4 (right).

Figure 5 .
Figure 5. Interval Method in an inclined borehole.Left: scheme of the test execution.Right: Shear wave velocity profile

Figure 7 .
Figure 7. Refracted ray path method based on Snell's law in an inclined borehole

Figure 9 .
Figure 9. Iterative process of the RRMI method.a) Number of iterations versus error.b) 2D ray-tracing layer model obtained for the 10th receiver in a 60° inclined borehole

Figure 10 .
Figure 10.Two-layer and three-layer synthetic soil model acquisition scheme

Figure 11 .Figure 12 .Figure 13 .
Figure 11.Shear wave velocity profiles obtained by DMI (corrected time plot), IMI, MIMI and RRMI methods along the well with a dip angle of 30° for two-layer models

Figure 14 .Figure 15 .Figure 16 .
Figure 14.Shear wave velocity profiles obtained by DMI (corrected time plot), IMI, MIMI and RRMI methods along the well with a dip angle of 30° for three-layer models

Figure 17 .Figure 18 .
Figure 17.Seismic velocity calculation error as a function of well dip angle

Figure 20 .
Figure 20.Shear wave velocity profiles used to model wave propagation.a) Shear wave velocity profiles from 2D MASW tests; the magenta boxes are the extracted sections for synthetic models.b) 2D synthetic models, the locations of the seismic source and receivers are marked with a red asterisk and black circles, respectively, along the well (dip angle 37°, length 30 meters).

Figure 21 .
Figure 21.Shear wave velocity profiles were calculated using the DMI (corrected time plot), IMI, MIMI, and RRMI methods.Top; synthetic records generated by wave propagation using FDM.Middle; corrected time (Tc) plot.Bottom; Shear wave velocity profiles along wells with 37° inclination for three models: M-01, M-02, and M-03.
).A 0.1 m grid was used with a source located 3.0 m from the wellhead.The borehole was 30 m long with a dip angle of 37°, and the receivers were spaced 1.0 m apart along the length of the borehole.This model generated synthetic wave propagation, and logs were obtained at each receiver, as shown in Figures 22-b and 22-c, respectively.

Figure 22 .
Figure 22. 2D synthetic model with complex soil layers used to obtain a shear wave velocity profile along a 37° dip well using the proposed methods.a) 2D synthetic model, the seismic source and receiver locations are indicated by red asterisks and green dots, respectively.b) Wave propagation at 4ms. c) Synthetic shot recorded by receivers along the well.d) 1D S-wave velocity profile calculated using the DMI, IMI, MIMI and RRMI methods.
.B. and Z.A.; methodology, P.B.; validation, P.B. and Z.A.; formal analysis, P.B.; writingoriginal draft preparation: P.B.; writing-review and editing: Z.A.All the authors have read and agreed to the published version of the manuscript.

Table 2 . Shear wave velocity and travel times were calculated in three-layer models in wells with dip angles of 30°, 45°, and 60°
Vs(m/s) Time (ms) Vs (m/s) Time (ms) Vs (m/s) Time (ms) Vs (m/s) Time (ms)