In viscoacoustic medium, the pressure p (x, t) and particle velocity v (x, t) satisfy a constitutive relation as ∂ p x , t ∂ t = ∫ − ∞ ∞ ∂ κ x , t − t ′ ∂ t ∇ ⋅ v x , t ′ d t ′ , (1) where v = (v _x , v _y , v _z ), and subscripts x, y and z denote particle velocity components along different axises. κ(x, t) is a time-dependent bulk modulus and can be constructed by superposing N groups of relaxation mechanisms (Fichtner and van Driel, 2014): κ x , t = κ 0 x 1 + 1 Q x ∑ l = 1 N D l ⁡ exp − t τ l H t , (2) where κ 0 ( x ) = ρ ( x ) v p 2 ( x ) is the relaxed bulk modulus, ρ(x) and v _p (x) denote the density and P-wave velocity, Q(x) is the quality factor, and H(t) is the step function. D ^(l) and τ ^(l) (1 < = l < = N) are the weights and decay times of different relaxation mechanisms, which can be computed by solving a non-linear optimization problem using a simulated annealing process to fit the frequency-independent Q in a limited band. In this study, we set the reference frequency to 1 Hz and the frequency band to [1, 150] Hz with three relaxation mechanisms. The resulting coefficients are shown in Table 1.

Taking the time derivative of κ(x, t), we have ∂ κ x , t ∂ t = κ 0 x 1 Q x ∑ l = 1 N − D l τ l exp − t τ l H t + 1 + s Q x δ t , (3) where δ(t) is the Kronecker delta function, and s is the summation of weights, i.e., s = ∑ l = 1 N D ( l ) . Inserting Eq. 3 into the constitutive relation yields ∂ p x , t ∂ t = κ 0 x 1 + s Q x ∇ ⋅ v x , t − 1 Q x ∑ l = 1 N Φ l x , t , (4) where the memory variables Φ^(l)(x, t) are defined as Φ l x , t = D l τ l ∫ − ∞ ∞ ⁡ exp − t − t ′ τ l H t − t ′ ∇ ⋅ v x , t ′ d t ′ , (5) and satisfy ∂ Φ l x , t ∂ t = − 1 τ l Φ l x , t + D l τ l ∇ ⋅ v x , t . (6)

Considering the second Newton law and assembling Eqs 4, 6, we obtain a viscoacoustic wave equation with the explicitly expressed Q as ρ x ∂ v x , t ∂ t = ∇ p x , t , ∂ p x , t ∂ t = ρ x v p 2 x 1 + s Q x ∇ ⋅ v x , t − 1 Q x ∑ l = 1 N Φ l x , t + f t δ x − x s , ∂ Φ l x , t ∂ t = − 1 τ l Φ l x , t + D l τ l ∇ ⋅ v x , t , (7) where x _s denotes the source location, and f(t) is the source time function.

Compared with existing viscoacoustic wave equations, such as the standard linear solid model (Carcione, 1990; Guo and Mcmechan, 2017) and fractional-Laplacian method (Zhu and Harris, 2014; Xing and Zhu, 2019), Eq. 7 has the following advantages in seismic modeling. 1) The quality factor Q is explicitly incorporated in the wave equation, which does not need to be converted to stress and strain relaxation times as in the standard linear solid method. 2) Eq. 7 can be efficiently solved using any time-domain finite-difference schemes, which does involves the Fourier transform or more complicated mixed-domain solvers as in the fractional Laplacian based methods. 3) As shown in Fichtner and van Driel (2014), the frequency-dependent Q effect can also be simulated by recomputing the weights D ^(l) and relaxation times τ ^(l).

To conform to the topographic surface of the Loess Plateau, here we use a vertically deformed grid to map the irregular physical space x = (x, y, z) to a regular computational domain α = (α, β, γ), which is given by α = x , β = y , γ = γ m a x z + ζ x , y z m a x + ζ x , y , (8) where ζ(x, y) is the elevation of topography (Figure 1), z _max is the maximum depth of the region of interest, γ _max = z _max +ζ _max , and ζ _max is the maximum of elevation.

With the coordinate transform relation in Eq. 8, the partial derivatives can be computed as ∂ ∂ x = ∂ ∂ α + γ m a x − γ z m a x + ζ x , y ∂ ζ x , y ∂ x ∂ ∂ γ , ∂ ∂ y = ∂ ∂ β + γ m a x − γ z m a x + ζ x , y ∂ ζ x , y ∂ y ∂ ∂ γ , ∂ ∂ z = γ m a x z m a x + ζ x , y ∂ ∂ γ . (9)

By setting c x = γ m a x − γ z m a x + ζ x , y ∂ ζ x , y ∂ x , c y = γ m a x − γ z m a x + ζ x , y ∂ ζ x , y ∂ y , c z = γ m a x z m a x + ζ x , y , (10)

The viscoacoustic wave equation in Eq. 7 can be rewritten in the new coordinate as ρ α ∂ v α , t ∂ t = ∇ ̄ p α , t , ∂ p α , t ∂ t = ρ α v p 2 α 1 + s Q α ∇ ̄ ⋅ v α , t − 1 Q α ∑ l = 1 N Φ l α , t + f t δ α − α s , ∂ Φ l α , t ∂ t = − 1 τ l Φ l α , t + D l τ l ∇ ̄ ⋅ v α , t , (11)

With a generalized partial derivative operator ∇ ̄ = ∂ ∂ α + c x ∂ ∂ γ ∂ ∂ β + c y ∂ ∂ γ c z ∂ ∂ γ . (12)

The proposed viscoacoustic wavefield simulation in a vertical stretched grid can be summarized into the following steps. First, with the elevation of topography ζ(x, y), we compute its derivative with respect to x and y and the topography related coordinate stretching parameters c _x and c _y . Then, the irregular physical space is mapped to a regular computational domain by vertical stretching according to Eq. 8, in which velocity, density and Q parameters onto the computational grids are calculated using a linear interpolation method. Next, the viscoacoustic wave equation 11 is solved in a regularly-sampled computational domain using any available numerical solvers. Finally, the wavefields in the real physical space are reconstructed by mapping the extrapolated wavefields back with the same interpolation algorithm used in the second step.

Considering the trade-off of computational accuracy and efficiency, we choose the time-domain staggered-grid finite-difference scheme to solve the viscoacoustic wave Eq. 11. Because of vertical stretching, it is difficult for the standard staggered-grid finite-difference method (Virieux, 1986) to accurately compute the partial derivatives such as ∂v _x /∂γ and ∂v _y /∂γ. Here we discretize and solve the wave equation using a fully staggered-grid approach (Lebedev, 1964; Lisitsa and Vishnevskiy, 2010; de la Puente et al., 2014). The nodes distribution is presented in Figure 2. Unlike the standard staggered-grid approach, the fully staggered-grid finite-difference scheme has three particle velocity components at the half grids (yellow triangles in Figure 2B), and three additional pressure wavefields (red circles in Figure 2B). For the time derivative, we still adopt a second-order finite-difference method, which results in the following updating scheme: p n + 1 = p n + Δ t ρ v p 2 1 + s Q D α + c x D γ v x n + 1 / 2 + D β + c y D γ v y n + 1 / 2 + c z D γ v z n + 1 / 2 − 1 2 Q ∑ l = 1 N Φ l , n + Φ l , n + 1 , Φ l , n + 1 = 1 − Δ t / 2 τ l 1 + Δ t / 2 τ l Φ l , n + Δ t 1 + Δ t / 2 τ l D l τ l D α + c x D γ v x n + 1 / 2 + D β + c y D γ v y n + 1 / 2 + c z D γ v z n + 1 / 2 , v x n + 3 / 2 = v x n + 1 / 2 + Δ t ρ D α p n + 1 + c x D γ p n + 1 , v y n + 3 / 2 = v y n + 1 / 2 + Δ t ρ D β p n + 1 + c y D γ p n + 1 , v z n + 3 / 2 = v z n + 1 / 2 + Δ t ρ c z D γ p n + 1 , (13) where D _i (i = α, β, γ) denotes the finite-difference operator for first-order spatial derivative. All components of pressure, particle velocity and memory variable wavefields are iteratively computed based on Eq. 13, and each partial derivative can be accurately calculated. Taking p ⁿ⁺¹ at (i, j, k) as an example, we compute D α v x n + 1 / 2 and D γ v x n + 1 / 2 using x-component particle velocities at (i ± m/2, j, k) and (i, j, k ± m/2), compute D β v y n + 1 / 2 and D γ v y n + 1 / 2 using y-component particle velocities at (i, j ± m/2, k) and (i, j, k ± m/2), and compute D γ v z n + 1 / 2 using z-component particle velocity at (i, j, k ± m/2), respectively. m = (1, 2 … M), where M denotes the half length of finite-difference operator. The other variables at different nodes can be calculated using the same finite-difference scheme. In numerical examples, we set the M = 4 and achieve an eighth-order accuracy in space, and the staggered-grid finite-difference coefficients are calculated using the Taylor expansion. To ensure the stability of the finite-difference solver, the temporal increment needs to satisfy the Courant-Friedrichs-Lewy (CFL) condition. In the vertically stretched grid, the relation between time and space increments is (the detailed derivation is given in Appendix A) Δ t ≤ Δ h c v m a x 1 + s Q m i n 1 + c x , m a x 2 + 1 + c y , m a x 2 + c z , m a x 2 , (14) where c is the summation of finite-difference coefficients, v _max denotes the maximum velocity, Q _min denotes the minimum quality factor value, and c x i , m a x (x _i = x, y, z) denotes the maximum value of topography associated coefficients in Eq. 10. For acoustic medium and a flat surface, Q tends to infinity, c _x = c _y = 0 and c _z = 1. Eq. 14 is simplified to a traditional CFL condition: Δ t ≤ Δ h 3 c v m a x . (15)

The free-surface boundary condition in the viscoacoustic modeling is implemented using a simple image method (Zhang et al., 2012). In addition, because two types of grids at different locations are coupled in the fully staggered-grid finite-difference scheme, spurious waves may occur if loading source incorrectly (Koene et al., 2021). To avoid these spurious waves, we adopt the strategy proposed by Lisitsa and Vishnevskiy (2010) and de la Puente et al. (2014), which is implemented for loading a point source by adding the source function at a node of (i, j, k) and at additional sub-nodes of (i ± 1/2, j ± 1/2, k) (i ± 1/2, j, k ± 1/2), (i, j ± 1/2, k ± 1/2) with a scale of 0.25. Seismograms are extracted using a similar strategy by summing the weighted records at a main node and sub-nodes of receiver locations.

A simulation example for a homogeneous model is presented in Figure 3. P-wave velocity is 3 km/s, Q is 50, and a 200-m-thick vacuum layer is set at the top to simulate surface reflections. Acoustic and two kinds of viscoacoustic solvers are used to compute seismic wavefields, and the result of generalized standard linear solid (GSLS) method is used as a reference for comparison. Because of the phase dispersion and energy dissipation, viscoacoustic waves propagate faster than acoustic waves and have much weaker amplitudes (Figures 3C,D). The surface at the z = 0 km produces a reflection with opposite polarity compared with the direct wave (Figure 3B). The waveforms computed using the proposed method have good agreements with those of GSLS in terms of both phases and amplitudes (red dot and black solid lines Figures 3C,D).

The first example is a three layer model, of which the velocity and Q values are shown in Figure 4A. A sinusoidal topography is designed to simulate the irregular surface. This model is discretized on a 601 × 501 grid with a 10-m spatial increment. A Ricker wavelet with the peak frequency of 15 Hz is used as the source time function. The source is deployed 10 m below the topography, and 501 receivers are evenly distributed horizontally on the observation surface. Time increment is 1 ms and duration is 6 s. Figure 4B shows the vertically distorted model in the computational domain. The topography has been fattened, but the originally horizontal interfaces in subsurface are curved. This coordinate conversion make it easy to simulate wave propagation in the near surface region.

The snapshots computed by solving the viscoacoustic and acoustic wave equations with and without the free surface are presented in Figures 5, 6. In the computational domain, the wavefronts of direct and reflected waves are distorted because of c _x , c _y, and c _z coefficients in Eq. 11 associated with the topography (Figure 5A–C, Figure 6A–C). In contrast, in the physical domain, the direct wave becomes regular with a half-circle wavefront, and transmitted and reflected waves are generated at the interfaces (Figure 5D–F, Figure 6D–F). The free surface at topography produces complicated boundary reflections, which intersects with deep effective wavefields (Figures 6D–F). Compared with acoustic modeling results, viscoacoustic wavefields have similar amplitudes and phases during early time, but they show considerably traveltime and amplitude differences as the propagation time increases (Figure 7).

The common-shot gathers are presented in Figures 8, 9. For comparison, we also compute the record using the GSLS method and utilize it as a reference, in which three relaxation mechanisms are used to approximate a frequency-independent Q. Without the free surface, common-shot records only have a direct wave and two subsurface reflections. The direct waves computed using acoustic and viscoacoustic simulations have similar waveforms due to short propagation time. By contrast, the deep reflectors of viscoacoustic records, especially for the one from second interface, have different phases and much weaker amplitudes than acoustic records (Figure 8D). The GSLS benchmark result has a good agreement with that of the proposed method, indicating that the new method can accurately simulate viscoacoustic propagation with irregular topography (black solid lines and red dot lines in Figures 8D, 9D). In addition, the reflected events are no longer hyperbolic due to topographic surface, and incorporating the free surface in the modeling introduces many surface-related multiples (Figure 9).

According to the outcrop of Chinese Loess Plateau, we build a realistic model for seismic modeling (Figure 10). Four layers are designed to mimic real loess structures: 1) a weathering layer with dry loess, 2) a low-velocity layer with wet loess, 3) a transition layer with water saturated loess, and 4) a Tertiary soil layer. Detailed P-wave velocity and Q values of these layers are given in Table 2. The maximum thickness of loess is about 300 m (Figure 10A), and low Q values in the dry and wet loess layers can produce strong seismic attenuation. This model is discretized onto a 1,000 × 2601 grid with a 5-m spacing. The source is located at x = 6.15 km horizontally and z = −1.18 km vertically with a 15-Hz Ricker wavelet as the source time function. The projected velocity in the computational domain is shown in Figure 10B, of which the surface is flattened and subsurface structures are slightly distorted.

The snapshots in the computational domain and physical space computed with and without the free surface are shown in Figures 11, 12. With a top absorption boundary condition, the reflected waves of subsurface interfaces can be clearly distinguished from direct waves (Figure 11). The wavefronts in the computational domain are slightly distorted and appear to be non-continuous because of irregular topography (left column in Figure 11). After mapping back to the physical domain, the wavefields, especially for reflections, have regular half-circle wavefronts (right column in Figure 11). By setting the topography as a free surface, seismic waves propagate back and forth between top surface and loess bottom vertically and between hill and ditch horizontally, resulting in strong multiple-like scatterings. These scatterings contaminate effective reflections generated from deep interfaces (Figure 12). Comparisons between acoustic and viscoacoustic modeling show a similar result as in the example of layer model (Figure 13). Their waveforms are similar when the propagation time is less than 0.5 s (Figure 13A). But as the time increases to 1 s, the difference of phases and amplitudes can be visually observed (Figure 13B). With a further longer duration, viscoacoustic wavefields have considerably weaker amplitudes than acoustic modeling results (Figures 13C,D).

The corresponding common-shot gathers are presented in Figure 14. Large velocity contrasts between loess layers and basement produce complicated refractions around the first break, and irregular topography results in distorted non-hyperbolic reflections (Figure 14A). Incorporating the free surface into viscoacoustic modeling produces strong scattering noises near the direct wave, and deep effective reflections are totally submerged in the noise (Figure 14B). For comparison, we also compute acoustic records using the same model and observation settings. Similar to the analysis for snapshots, the events at a short time have similar phases and amplitudes on acoustic and viscoacoustic records (Figures 14C,D). But considerable difference of deep reflections can be observed as the time increases due to attenuation-related phase dispersion and amplitude dissipation. In addition, large-offset refractions and scatterings near direct wave of viscoacoustic record are much weaker than those of acoustic modeling results. This is caused by strong attenuation in the low-Q loess layers.

Because 2D modeling cannot accurately simulate the amplitudes of seismic waves, we apply the proposed viscoacoustic modeling method to a more realistic 3D model (Figure 15). Compared with 2D Loess Plateau model, geographic characteristics of loess distribution, including typical ditch, hill and valley, are more clearly visible in the 3D model. The setting for Q in shallow loess layers is the same to that in 2D model. For deep sandstone strata, the well-logging data and rock physics analysis show that P-wave velocity and Q models crudely satisfy an empirical relation as Q = v p , where v _p denotes P-wave velocity with the unit of m/s. This model is discretized by 475 × 1,301×800 nodes with a 10-m increment. 650 × 200 receivers along inline and crossline directions are uniformly deployed to record seismograms (blue dots in Figure 16). Three experiments are carried out with sources in the ditch, hill and valley, respectively (red stars in Figure 16). Time increment in the finite-difference is 0.5 ms and duration is 7 s.

The snapshots and common-shot gathers are presented in Figures 17, 18. Similar to 2D modeling, the interactions between topographic free-surface and loess bottom produces strong scatterings in near-surface region, which contaminate deep reflections (Figure 17). This leads to typical “black-triangle” noise below the source location in the records (Figure 18A). Without subsequent processing, it is difficult to visually find hyperbolic reflections. Because of the geometric spreading in a 3D volume, wavefield amplitude decays faster during propagation than in the 2D case. In addition, seismic attenuation introduces strong amplitude dissipation for viscoacoustic waves (right column in Figure 17). This results in rather weak amplitudes at large-offset traces and large observation time compared with acoustic recordings (Figures 18B,C).

The modeling results of the other two experiments are shown in Figures 19–21. Thicker loess sediments in the hill and valley produce stronger attenuation effect on seismic waves than that in the ditch, resulting in large amplitude and phase difference between acoustic and viscoacoustic snapshots (Figure 19). The phase dispersion and amplitude dissipation of viscoacoustic recordings become more pronounced. The waveforms with large offsets and times, including near-surface-related scatterings and deep reflections, are attenuated quickly (Figures 20, 21). In contrast to acoustic recordings, no events can be observed visually at the time greater than 5 s. The low SNR and weak amplitudes of recordings pose a great challenge for subsequent seismic data processing and imaging.

Figure 22 shows the snapshots of 1.5 s on the observation surface for three experiments. In the first experiment, strong scatterings are generated and propagate along the north-south ditch (Figures 22A,B). This is because the loess in the ditch was eroded, and the outcrop has a relative higher velocity than the loess. In contrast, the experiments excited on the hill and valley show an isotropic scattering pattern (Figures 22C–F). In addition, the wavefields excited by the source on the hill are affected by attenuation more largely than those excited in the ditch and valley (Figures 22D,F). Comparisons of acoustic modeling result, viscoacoustic modeling result and a common-shot gather acquired from Chinese Loess Plateau are presented in Figure 23. Note that both acoustic and viscoacoustic modeling can simulate strong black-triangle scattering noise below the direct waves, which has a consistent kinematic pattern with the field data. But acoustic records have relatively stronger amplitudes for the time greater than 3 s (Figure 23B). By incorporating seismic attenuation, viscoacoustic modeling produces weak amplitudes for large-duration and large-offset events, which have a good agreement with the amplitude characteristic of field data (Figure 23C).