Viscoelastic Model and Synthetic Seismic Data of Eastern Rub’Al-Khali

Abstract

The Rub’ Al-Khali basin in Saudi Arabia remains unexplored and lacks data availability due to its remoteness and the challenging nature of its terrain. Thus far, there are neither digital geologic models nor synthetic seismic data from this specific area accessible for testing research techniques and analysis. In this study, we build a 2D viscoelastic model of the eastern part of the Rub’ Al-Khali basin and generate a corresponding dual-component seismic data set. We compile high-resolution depth models of compressional- and shear-wave velocities, density, as well as compressional- and shear-wave quality factors from published data. The compiled models span Neoproterozoic basement up to Quaternary sand dunes. We then use the finite-difference technique to model the propagation of seismic waves in the compiled viscoelastic medium of eastern Rub’ Al-Khali desert. In particular, we generate vertical and horizontal components of the shot gathers with accuracy to the fourth and second orders in space and time, respectively. The viscoelastic models and synthetic seismic datasets are made available in an open-source site for prospective re-searchers who desire to use them for their research. Users of these datasets are urged to make their findings also accessible to the geoscience community as a way of keeping track of developments related to the Rub’ Al-Khali desert.

1. Introduction

Numerical modeling of seismic wave propagation is a powerful way to understand real seismic data acquired in an area, because these methods closely simulate seismic waves from their generation at the source until recording at the receiver. Seismic numerical modeling is best done in a 3D sense in order to closely simulate wave propagation in 3D earth models. However, it can also be done in 1D, 2D, and 2.5D models in some applications to reduce the simulation time. Various methods are used to simulate 3D seismic wave propagation including finite-difference, e.g., [1,2,3,4], finite-volume, e.g., [5,6,7], and Lattice–Boltzmann, e.g., [8,9,10]. In particular, the finite-difference method covers an intuitive and practical modeling approach because it can be used for single-and multi-dimensional spaces and wave equations of any complexity. In the literature, the viscoelastic models and synthetic seismic data, calculated using finite-difference method, have been used as input in various studies [11,12,13,14]. Regardless of the simulation method, increasing levels of wave equation complexity can be used to model seismic wave propagation depending on how realistic the depth model is. The simplest but least realistic model uses the acoustic wave equation where only P-wave and density are considered. On the other hand, the most complex and realistic model uses an anisotropic poro-viscoelastic wave equation, which addresses almost all types of seismic waves and phenomena [15].

The highest level of modeling is not practical because of the need to know many generally unknown parameters (i.e., anisotropy) about subsurface layers, in addition to being very computationally intensive. The viscoelastic wave equation offers a practical compromise, especially in large subsurface models spanning a wide range of velocities. In addition to modeling P-waves, S-waves, and surface waves, viscoelasticity also addresses an elastic wave behavior. The Rub’ Al-Khali basin is situated below the Quaternary sand dunes and occupies most of the southern part of Saudi Arabia. It covers a large area that extends approximately 1300 km in the east–west direction and 300 km in the north–south direction extending from southwest Saudi Arabia toward the United Arab Emirates (UAE) [16] (Figure 1A). This basin is divided by the Qatar arch into two sub-basins: East and West Rub’ Al-Khali sub-basins (Figure 1B). In comparison with the UAE part, Rub’ Al-Khali basin in Saudi Arabia remains unexplored and lacks openly available data due to its remoteness and the challenging desert conditions. Most of the hydrocarbon fields in Saudi Arabia in the Rub’ Al-Khali basin are mainly located on its northeastern portion, including Shaybah, Ramlah, North, and South Kidan fields [17]. Alfara et al. [18] described the main challenges and uncertainties in the interpretation of 3D seismic data in the Shaybah field due to the complex and rugged topography. They successfully improved the quality of the data and resolved the uncertainties by combining 3D seismic data and borehole seismic measurements including vertical seismic pro ling (VSP) surveys. Stewart et al. [19,20] carried out a study on the Mesozoic petroleum system of the Rub’ Al-Khali basin, particularly the Triassic to Jurassic siliciclastic reservoirs. They also described the structural evolution of the Rub’ Al-Khali basin from the Precambrian to Neogene using a combination of regional seismic, gravity, and magnetic data.

Thus far, there are no publicly-available synthetic seismic data nor digital geological models of the Rub’ Al-Khali basin. In this study, we build a 2D viscoelastic model of the eastern Rub’ Al-Khali, which is close to the border of the UAE (Figure 1B) and generate the associated dual-component seismic data sets. Similar models and synthetic seismic data sets from the northwest, eastern, and central parts of Saudi Arabia have been made publicly available by the authors [21,22,23]. Models of P-wave velocity (V_p), S-wave velocity (V_s), Density (ρ), P-wave quality factor (Q_p), and S-wave quality factor (Q_s) that cover Neoproterozoic basement to Quaternary sand are first compiled from available literature [24]. These models are subsequently utilized in generating the 2D synthetic viscoelastic dual-component seismic data sets that include vertical component data and horizontal component data recorded by surface geophones using a conventional 2D acquisition geometry. The horizons, models, and synthetic seismic data sets are made publicly available at the following link: (https://www.dropbox.com/sh/9t099hhn1h8nea8/AAChpqVg3R_km2MlEMJ0SuxLa?dl=0).

The synthetic seismic data sets are saved in SEGY files, while horizons and faults are saved in a CSV text file.

2. Geology of the Rub’ Al-Khali Basin

Due to lack of outcrop exposure (except in the west margin), the Rub’ Al-Khali basin evolution is known only from geophysical data. Stewart [20] reported the process of development of structures in the Rub’ Al-Khali basin from the Late Precambrian to Neogene observed on seismic, magnetic, gravity, and well data. This data shows a relatively complete Phanerozoic succession with north–south and northwest–southeast trending basement structure (Figure 1B) that are formed during East African Antarctic Orogeny. The structural trend in the Rub’ Al-Khali basin has been affected and reactivated during several periods including Ediacaran—Cambrian Ara/Hormuz salt, Devonian to Late Carboniferous Hercynian orogeny, early Mesozoic Gondwana break up, and finally in the Late Cretaceous and Tertiary time due to effects of the Arabian plate colliding with the Eurasian plate, resulting in the emplacement of ophiolites in the Oman Mountains.

3. Development of Digital Depth Models

There are almost no published data about the eastern parts of the Rub’ Al-Khali basin in Saudi Arabia. Only limited published data on the geological or geophysical properties were available for the Rub’ Al-Khali basin in Saudi Arabia. In this study, we performed an extensive literature review to compile one geophysical and geological model in the eastern Rub’ Al-Khali area using references specific to Rub’ Al-Khali. In cases where data are not available in these references, we obtain the missing data from references studying the same formations in nearby areas of Saudi Arabia, UAE, or Oman. In terms of geologic column in the Rub’ Al-Khali basin, we found the following twenty-three formations: Quaternary (F-1), Hofuf Dam Hadrukh (F-2), Dammam (F-3), Rus (F-4), Umm Er Radhuma (F-5), Aruma (F-6), Wasia (F-7), Shuaiba (F-8), Biyadh (F-9), Hith (F-10), Arab (F-11), Hanifa & Tuwaiq Mountain (F-12), Dhruma (F-13), Marrat (F-14), Minjur (F-15), Jilh (F-16), Sudair (F-17), Khuff (F-18), Unayzah (F-19), Qusaiba (F-20), Qasim (F-21), Saq (F-22), and Basement (F-23) (Figure 2). Eventually, twenty-three layers and two major faults were defined and digitized at variable intervals (Figure 2). The Phanerozoic stratigraphy column of the eastern Arabian Peninsula is given in Figure 3 [25]. The stratigraphy of layers and faults of the model were built mostly based on the geologic information provided by [19]. The elastic properties such as densities and velocities were determined using median values from digitized well-log data provided by several published articles as reported in

Table 1. Lithology, densities, compressional and shear-wave speeds, compressional and shear wave quality factors of layers of the eastern Rub’ Al-Khali basin model.

Layer	Formation	Lithology	ρ (kg/m3)	Vp (m/s)	Vs (m/s)	Qp	Qs
F-1	Quaternary	Sand	1500 a	850 a	501 a	29.15 i	22.36 i
F-2	Hofuf Dam Hadrukh	Mixed limestone and siliciclastic	1877 b	1835 b	1180 c	42.84 i	34.35 i
F-3	Dammam	Limestone	2289 b	3110 b	1870 d	55.77 i	43.24 i
F-4	Rus	Anhydrite	2397 b	5252 b	2470 d	65.32 i	49.70 i
F-5	Umm Er Radhuma	Limestone	2028 b	3318 b	1978 d	57.60 i	44.47 i
F-6	Aruma	Limestone	2037 b	3558 b	1672 d	52.25 i	40.89 i
F-7	Wasia	Sandstone	2277 b	3233 b	2142 c	56.86 i	46.28 i
F-8	Shuaiba	Limestone	2037 b	3010 b	1818 d	54.86 i	42.63 i
F-9	Biyadh	Sandstone	2364 b	4045 b	2701 c	63.60 i	51.97 i
F-10	Hith	Anhydrite	2874 e	4483 e	2328 e	66.96 i	48.24 i
F-11	Arab	Limestone	2400 e	5399 e	2748 e	73.47 i	52.42 i
F-12	Hanifa & Tuwaiq Mountain	Limestone	2550	5698 e,k	2903 e	75.48 i	53.88 i
F-13	Dhruma	Limestone	2458 f,k	5033 f,k	2870 d	70.94 i	53.57 i
F-14	Marrat	Shale	2410 f	3272 f	1436 c	57.20 i	37.89 i
F-15	Minjur	Sandstone	2394 f	3930 f	2499 c	62.69 i	49.99 i
F-16	Jilh	Dolomite	2400 f,k	4823 f,k	2761 d	69.45 i	52.54 i
F-17	Sudair	Shale	2372 g,k	5182 g,k	2674 g	71.99 i	51.71 i
F-18	Khuff	Dolomite	2639 g,k	4953 g,k	2530 g	70.38 i	50.3 i
F-19	Unayzah	Sandstone	2405 g	3752 g	2085 g	61.25 i	45.66 i
F-20	Qusaiba	Shale	2486 g	3898 g	2143 g	62.43 i	46.29 i
F-21	Qasim	Sandstone	2380 h	3685 h	2453 c	60.70 i	49.53 i
F-22	Saq	Sandstone	2350 h	3765 h	2508 c	61.36 i	50.08 i
F-23	Basement	Igneous and metamorphic	2800 j	6380 j	3580 j	79.87 i	59.83 i

^a [28], ^b [18], ^c [22,23], ^d [26], ^e [29], ^f [30], ^g [31], ^h [32], ⁱ [27], ^j [33], ^k [34].

[18,22,23,26,27,28,29,30,31,32,33,34], and available Vp and density data were plotted in Figure 4. If logs are not available, values for these parameters from nearby areas are used. Due to the lack of availability of S-wave velocities for some formations, we use Vp-Vs relations to complete the whole stratigraphic column. We used the relation developed by [26] for carbonate lithologies, and relations presented by [22] for sandstones and shales. The remaining properties, which are P- and S-wave quality factors (Qp and Qs), were calculated using a method proposed by [27], which involves taking the square root of the corresponding value of P- and S-wave, respectively. The complete geophysical properties including P-wave velocity, S-wave velocity, density, lithology, P- and S- wave quality factors of each formation are listed in Table 1. Moreover, Figure 5 shows plots of the compiled Vp, Vs, ρ, Qp, and Qs models.

4. The Synthetic Seismic Data Generation Procedure

We apply 2D finite-difference approximation methodology to calculate the synthetic seismic data, using the available viscoelastic model inputs. This method involves the application of 2D wave equation to simulate the propagation of seismic waves in viscoelastic media that is determined by a modified Newton’s equation of motion and Hooke’s law, as described in [35,36]. We propagate the waves in a 2D viscoelastic medium by computing the 3 strain tensor components along the x-axis and z-axis, and also compute their associated stress tensors. Then, we calculate the divergence of the stresses which generate accelerations that are further extrapolated to produce displacements to succeeding time steps. The finite-difference method used here utilizes staggered grids to compute the spatial derivatives using the centered finite-difference approximations. The finite-difference algorithm inputs the viscoelastic properties (i.e., Vp, Vs, Qp, and Qs) from the compiled models to compute Lame’s coefficients and the wave’s relaxation parameters.

The finite-difference viscoelastic scheme begins with generation and propagation of the source at an initial time t = 0. Then it computes the wave-speeds using the finite-difference scheme by increasing time in little time steps while modifying the wave-speed in both x and z-directions, respectively. Subsequently, this process updates the stress field and makes it staggered in time relative to particle wave-speed fields. The process is repeated with another time step, updating the particle wave-speeds in this case until it covers the entire grid. At the boundaries of the grid, we utilize absorbing boundary conditions to minimize the effects of reflections. Usually, the computation of the finite-difference solutions relies on the frequency band of the source wavelet. In this case, we apply a typical Ricker wavelet having a zero-phase and a 20-Hz peak frequency. The finite-difference solution of the wave equation is usually associated with numerical dispersion if appropriate grid spacing is not used. To avoid this problem, we apply the general formula for preventing dispersion in the finite-difference scheme with accuracy close to spatial fourth-order and close to second order in time. This means that we must have a minimum of five grid points for every wavelength within a viscoelastic finite-difference scheme [37]:

$$\Delta \mathrm{x} = \Delta \mathrm{z} \le \frac{{\mathrm{V}}_{\min }}{5{\mathrm{f}}_{\max }}=1.5 \mathrm{m}$$

(1)

The grid sizes $\Delta$x (along x-axis) and $\Delta$z (along z-axis) are assumed to be equal, Vmin is the minimum wave speed, and fmax is the maximum frequency. In our model, Vmin = 501 m/s corresponding to the minimum shear-wave speed and fmax = 60 Hz, which is the frequency at which the wavelet amplitude spectrum drops below 1% of its maximum value. Consequently, we calculate a value of x = z = 1.5625 m, which is approximated to 1.5 m. Besides the dispersion criteria, the time-marching step of the finite-difference also needs to fulfill the Courant-Friedrichs-Lewy (CFL) stability criteria for convergence in 2D media [37], which is given as:

$$\Delta \mathrm{t} \le \Delta \mathrm{z} \frac{0.606\Delta \mathrm{x}}{{\mathrm{V}}_{\max }}=0.14 \mathrm{msec}$$

(2)

where Vmax is the maximum wave speed in the model. We use the basement compressional-wave speed as the maximum velocity (i.e., Vmax = 6380 m/s), then we compute the time step as t = 0.142853 ms, which is approximated to 0.14 ms. All these parameters are embedded in the 2D finite-difference wave field modelling open-source code [37], which we utilized for generating the viscoelastic synthetic seismic data. To satisfy the boundary conditions, we used absorbing criteria with a taper of 375 grid points at the bottom, right, and left edges of the model and a free-surface boundary condition for the top margin of the model. We compute the taper using the formula described in [37]:

$$\mathrm{Taper} = \frac{4{\mathrm{V}}_{\max }}{\Delta \mathrm{x}\times {\mathrm{f}}_{\max }}=375$$

(3)

Regarding the data acquisition, we apply a conventional 2D seismic land survey geometry. Starting from point x = 0 m, sources cover the whole horizontal length to point x = 40,000 m with a source interval of 50 m. Similarly, receivers cover a horizontal length from point x = 0 m to point x = 40,000 m, but with a receiver separation of 25 m. Receivers and shots are placed at 15 m depth below the top interface [37] to avoid the free-surface effect. The offset ranges from −40,000 to 40,000 m, which is desirable for most seismic data analysis and processing such as; velocity analysis, anisotropy, and amplitude variation with offset (AVO) analysis. Parameters used for the generation of synthetic seismic data sets are tabulated and presented in

Table 2. List of parameters used to construct synthetic seismic data.

Parameter	Value
Source wavelet	20 Hz zero-phase Ricker
Time sampling interval for finite-difference calculation	0.14 ms
Square-grid size for finite-difference calculation	1.5 m
Receiver spacing	25 m
Shot spacing	50 m
Recording time sampling interval	4.0 ms
Total recording time	6.0 s
Total number of receivers	1535
Receivers x-axis	0 to 40,000 m
Receivers z-axis	15 m
Shots x-axis	0 to 40,000 m
Shots z-axis	15 m
Offset range	−40,000 to 40,000 m
Total number of shots	801 shots

. We generated both horizontal and vertical components of the synthetic viscoelastic seismic data sets. Samples of the vertical component and horizontal component of displacements recorded at the beginning, center, and end of the model are presented in Figure 6 and Figure 7, respectively.

For generating the synthetic seismic data, we used a PC workstation running a Linux operating system with an Intel^® Xeon^® Processor E5-2687 3.0 GHz 10C and a memory of 256 GB 2133 MHz DDR4 ECC RDIMM. A single shot gather consumed approximately 36 h to generate on this workstation. We attribute this length to the relatively low minimum velocity and high maximum frequency of our model.

5. Conclusions

We have put together a 2D viscoelastic model from the eastern Rub’ Al-Khali basin and produced the associated dual-component synthetic seismic data sets. Regardless of the large number of layers involved in the model, the data appears to have good quality in capturing most important reservoirs including the relatively deep Khuff Formation. We provided the geologic model and associated seismic data sets, available online publicly for those who want to test their methods, and we encourage them to make their output open publicly too.