# Viscoelastic Model and Synthetic Seismic Data of Eastern Rub'Al-Khali
## 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 [PERSON]'s equation of motion and [PERSON]'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 [PERSON]'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  
Figure 5: (a) Vp model; (b) Vs model; (c) density model; (d) Qp model; (e) Qs model. 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}\ \leq\ \frac{\mathrm{V_{min}}}{5 \mathrm{f_{max}}}=1.5\ \mathrm{m} \tag{1}\]  
The grid sizes \(\Delta\mathrm{x}\) (along x-axis) and \(\Delta\mathrm{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}\ \leq\ \Delta\mathrm{z}\ \frac{0.606\Delta\mathrm{x}}{ \mathrm{V_{max}}}=0.14\ \mathrm{msec} \tag{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.