Viscoacoustic Wave Equation Research Articles

Modeling viscoacoustic wave propagation using a new spatial variable-order fractional Laplacian wave equation

We have developed a new time-domain viscoacoustic wave equation for simulating wave propagation in anelastic media. The new wave equation is derived by inserting the complex-valued phase velocity described by the Kjartansson attenuation model into the frequency-wavenumber domain acoustic wave equation. Our wave equation includes one second-order temporal derivative and two spatial variable-order fractional Laplacian operators. The two fractional Laplacian operators describe the phase dispersion and amplitude attenuation effects, respectively. To facilitate the numerical solution for our wave equation, we use the arbitrary-order Taylor series expansion (TSE) to approximate the mixed-domain fractional Laplacians and achieve the decoupling of the wavenumber and the fractional order. Then, our viscoacoustic wave equation can be directly solved using the pseudospectral method. We adopt a hybrid pseudospectral/finite-difference method (HPSFDM) to stably simulate wave propagation in arbitrarily complex media. We validate the high accuracy of our approximate dispersion term and approximate the dissipation term in comparison with the accurate dispersion term and accurate dissipation term. The accuracy of the numerical solutions is evaluated by comparison with the analytical solutions in homogeneous media. Theory analysis and simulation results indicate that our viscoacoustic wave equation has higher precision than the traditional fractional viscoacoustic wave equation in describing constant- Q attenuation. For a model with Q < 10, the calculation cost for solving the new wave equation with TSE HPSFDM is lower than that for solving the traditional fractional-order wave equation with TSE HPSFDM under the high numerical simulation precision. Furthermore, we examine the accuracy of HPSFDM in heterogeneous media using several numerical examples.

Q-compensated reverse-time migration for natural gas hydrate using fractional viscoacoustic wave equation

Natural gas hydrate (NGH) is a potential clean alternative energy source for fossil fuels. In seismic imaging profiles, NGH is often identified by the bottom-simulating reflection (BSR), which is characterized by strong reflection amplitude and negative polarity. High-resolution and amplitude-preserved seismic imaging are demanded for the detection of NGH. However, the traditional acoustic reverse-time migration ignores the attenuation characteristics of the medium, which leads to reduced amplitude and distorted phase of the seismic wave. In particular, when the hydrate saturation is low or the underlying formation does not contain free gas, it is difficult to observe the identifiable BSR in the traditional acoustic imaging profile, which causes difficulties in the identification of NGH. Here, we introduce the fractional viscoacoustic wave equation to perform the Q-compensated reverse-time migration (Q-RTM) for NGH, which can accurately recover the amplitude loss, correct the phase distortion, and provide high-resolution and high-illumination imaging results. Finally, Q-RTM can effectively enhance the BSR, reduce the uncertainty of hydrate identification, help to confirm the location and spatial distribution of the gas hydrate-bearing sediments, further refine the geological properties, and provide some theoretical basis for the exploitation and drilling of hydrate.

Constant-Q wave propagation and compensation by pseudo-spectral time-domain methods

Seismic attenuation is usually described by a constant-Q (CQ) model that assumes the seismic quality factor (Q) is independent of the frequency. To simulate the attenuation behaviors of seismic waves, we develop a pseudo-spectral time-domain (PSTD) method to solve a CQ viscoacoustic wave equation. This method is nearly fourth-order accurate in time. Compared to the conventional temporal second-order PSTD method, the new PSTD scheme is verified to be more efficient. In some applications such as Q-compensated reverse-time migration (Q-RTM) and time-reversal imaging, one requires to simulate an anti-attenuation process. To realize this purpose, we switch our viscoacoustic PSTD modeling scheme into an amplitude-compensated PSTD modeling scheme by flipping the signs of the operators that dominate the amplitude loss. To control the numerical instability caused by high-frequency overcompensation in the amplitude-boosted modeling, we integrate a time-variant filter to the PSTD modeling scheme to suppress the high-frequency noise. Wavefield simulation examples in homogeneous media verify the temporal accuracy of our nearly fourth-order PSTD modeling scheme. A Q-RTM test of synthetic data is also presented to demonstrate the robustness of our amplitude-compensated PSTD modeling scheme.

Visco-acoustic full waveform inversion using decoupled fractional Laplacian constant-Q wave equation and optimal transport-based misfit function

Full waveform inversion (FWI) has the potential to recover high-resolution physical parameters of the Earth from seismic data. Considering the attenuation effects of the subsurface, we develop a new visco-acoustic FWI method in the time domain with a known Q model. The seismic wavefield is modelled by solving a decoupled fractional Laplacian visco-acoustic wave equation, which can better describe amplitude attenuation and phase distortion during wave propagation. To weaken the initial model dependence of visco-acoustic FWI, we introduce the optimal transport-based misfit functional which can measure the data residual globally. Numerical examples on the 2D Marmousi model demonstrate the superiority of the proposed method compared with the L2 norm-based method which measures the misfit of observed and synthetic seismic data locally, generating cycle-skipping issue. The inversion results of seismic data without low-frequency components further demonstrate the validity and effectiveness of the proposed method, which is insensitive towards low-frequency components, and may achieve reasonable inversion results for field data application.

An analytic signal-based accurate time-domain viscoacoustic wave equation from the constant-Q theory

We have introduced a concise time-domain wave equation to accurately simulate wave propagation in viscoacoustic media. The central idea behind this work is to dismiss the negative frequency components from a time-domain signal by converting the signal to its analytic format. The negative frequency components of any analytic signal are always zero because we can construct the viscoacoustic wave equation to honor the relaxation property of the medium for positive frequencies only. The newly proposed complex-valued wave equation (CWE) represents the wavefield with its analytic signal, whose real part is the desired physical wavefield, whereas the imaginary part is the Hilbert transform of the real component. Specifically, this CWE is accurate for weakly and strongly attenuating media in terms of the dissipation and dispersion, and the attenuation is precisely linear with respect to the frequencies. The CWE easily and flexibly models dispersion-only, dissipation-only, or dispersion-plus-dissipation seismic waves. We have verified these CWEs by comparing the results with analytical solutions, and we achieved nearly perfect matches. When extending to heterogeneous media, the results are consistent with those computed from the nonstationary operator-based Fourier integral method and from the standard linear solid equations.

Fractional Time Derivative Seismic Wave Equation Modeling for Natural Gas Hydrate

Simulation of the seismic wave propagation in natural gas hydrate (NGH) is of great importance. To finely portray the propagation of seismic wave in NGH, attenuation properties of the earth’s medium which causes reduced amplitude and dispersion need to be considered. The traditional viscoacoustic wave equations described by integer-order derivatives can only nearly describe the seismic attenuation. Differently, the fractional time derivative seismic wave-equation, which was rigorously derived from the Kjartansson’s constant-Q model, could be used to accurately describe the attenuation behavior in realistic media. We propose a new fractional finite-difference method, which is more accurate and faster with the short memory length. Numerical experiments are performed to show the feasibility of the proposed simulation scheme for NGH, which will be useful for next stage of seismic imaging of NGH.

Source Wavefield Reconstruction in Fractional Laplacian Viscoacoustic Wave Equation-Based Full Waveform Inversion

We develop a fractional Laplacian viscoacoustic wave equation-based full waveform inversion (FWI) method. The main novelty is efficient reconstruction of the source wavefields in gradient computation by the adjoint state method. Our FWI is based on Fourier pseudo-spectral time-domain (PSTD) numerical solutions of the fractional Laplacian viscoacoustic wave equation, which can describe the frequency independent <inline-formula xmlns:mml="http://www.w3.org/1998/Math/MathML" xmlns:xlink="http://www.w3.org/1999/xlink"> <tex-math notation="LaTeX">${Q}$ </tex-math></inline-formula> (quality factor) behaviors of seismic waves accurately. The presented wavefield reconstruction strategy utilizes reverse time-marching formulae to compute the source wavefields backward in time, including an implicit formula in the interior domain and an explicit one in the perfectly matched layer (PML) absorbing domains. Since the wavefields in the entire domain are recovered, our method does not require storing massive boundary values like conventional reconstruction methods. To avoid numerical instability in reverse time reconstruction, we design a global-local checkpointing technique. When reverse time reconstruction encounters numerical instability, the forward propagation restarts from the nearest global checkpoint and proceeds to the unstable time step. Then, the reverse time reconstruction continues. The local checkpoints help to delay the numerical instability in the PML domains. Numerical examples verify the feasibility of our reconstruction strategy and the efficiency gain achieved by using this strategy in FWI.

Q-compensated reverse time migration in tilted transversely isotropic media

Anisotropy and absorption are critical to the modeling and analysis of seismic amplitude,phase, and traveltime data. To neglect any of these phenomena, which are often bothoperating simultaneously, degrades the resolution and interpretability of migrated images.However, a full accounting of anisotropy and anelasticity is computationally complex andexpensive. One strategy for accommodating these aspects of wave propagation, while keepingcost and complexity under control, is to do so within an acoustic approximation. Weset up a procedure for solving the time-domain viscoacoustic wave equation for tilted transverselyisotropic (TTI) media, based on a standard linear solid model and, from this, developa viscoacoustic reverse time migration (Q-RTM) algorithm. In this approach, amplitudecompensation occurs within the migration process through a manipulation of attenuationand phase dispersion terms in the time domain differential equations. Specifically, theback-propagation operator is constructed by reversing the sign only of the amplitude lossoperators, but not the dispersion-related operators, a step made possible by reformulatingthe absorptive TTI equations such that the loss and dispersion operators appear separately.The scheme is tested on synthetic examples to examine the capacity of viscoacoustic RTM to correct for attenuation, and the overall stability of the procedure.

Least-squares reverse time migration of multiples in viscoacoustic media

The migration of multiples can provide complementary information about the subsurface, but crosstalk artifacts caused by the interference between different-order multiples reduce its reliability. To mitigate the crosstalk artifacts, least-squares reverse time migration (LSRTM) of multiples is suggested by some researchers. Multiples are more affected by attenuation than primaries because of the longer travel path. To avoid incorrect waveform matching during the inversion, we propose to include viscosity in the LSRTM implementation. A method of LSRTM of multiples is introduced based on a viscoacoustic wave equation, which is derived from the generalized standard linear solid model. The merit of the proposed method is that it not only compensates for the amplitude loss and phase change, which cannot be achieved by traditional RTM and LSRTM of multiples, but it also provides more information about the subsurface with fewer crosstalk artifacts by using multiples compared with the viscoacoustic LSRTM of primaries. Tests on sensitivity to the errors in the velocity model, the Q model, and the separated multiples reveal that accurate models and input multiples are vital to the image quality. Numerical tests on synthetic models and real data demonstrate the advantages of our approach in improving the quality of the image in terms of amplitude balancing and signal-to-noise ratio.

Domain Decomposition for Large-Scale Viscoacoustic Wave Simulation Using Localized Pseudo-Spectral Method

Wave simulation in absorptive media using decoupled fractional Laplacian wave equation has received widespread attention in recent years, largely due to its precise description of frequency independent $Q$ and easy attenuation compensation in seismic processing. With many algorithms to solve the fractional Laplacian, $k$ -space pseudo-spectral method is predominantly used in academia, where the computing facilities support fast Fourier transforms. However, its global nature prevents the parallelization and computational efficiency of the forward solver, especially for large-scale applications. We propose to solve viscoacoustic wave equation using domain decomposition. A local Fourier basis is constructed around the truncated area to improve the periodicity and smoothness of the decomposed wave information. After independently simulating in the subdomains, a pointwise patching procedure is applied to maintain the continuity of the wavefield between subdomains. Numerical experiments show that this new algorithm obtains high computational efficiency without compromising the numerical stability condition of the traditional pseudo-spectral method. This work becomes more attractive for seismic inversion and imaging problems by improving its parallelization.

Arbitrary‐order Taylor series expansion‐based viscoacoustic wavefield simulation in 3D vertical transversely isotropic media

ABSTRACTTaking the anisotropy of velocity and attenuation into account, we investigate the wavefield simulation of viscoacoustic waves in 3D vertical transversely isotropic attenuating media. The viscoacoustic wave equations with the decoupled amplitude attenuation and phase dispersion are derived from the fractional Laplacian operator and using the acoustic approximation. With respect to the spatially variable fractional Laplacian operator in the formulation, we develop an effective algorithm to realize the viscoacoustic wavefield extrapolation by using the arbitrary‐order Taylor series expansion. Based on the approximation, the mixed‐domain fractional Laplacian operators are decoupled from the wavenumbers and fractional orders. Thus, the viscoacoustic wave propagation can be conveniently implemented by using a generalized pseudospectral method. In addition, we perform the accuracy and efficiency analyses among first‐, second‐ and third‐order Taylor series expansion pseudospectral methods with different quality factors. Considering both the accuracy and computational cost, the second‐order Taylor series expansion pseudospectral method can generally satisfy the requirements for most attenuating media. Numerical modelling examples not only illustrate that our decoupled viscoacoustic wave equations can effectively describe the attenuating property of the medium, but also demonstrate the accuracy and the high robustness of our proposed schemes.

Estimating P Wave Velocity and Attenuation Structures Using Full Waveform Inversion Based on a Time Domain Complex‐Valued Viscoacoustic Wave Equation: The Method

AbstractTo complement velocity distributions, seismic attenuation provides additional important information on fluid properties of hydrocarbon reservoirs in exploration seismology, as well as temperature distributions, partial melting, and water content within the crust and mantle in earthquake seismology. Full waveform inversion (FWI), as one of the state‐of‐the‐art seismic imaging techniques, can produce high‐resolution constraints for subsurface (an)elastic parameters by minimizing the difference between observed and predicted seismograms. Traditional waveform inversion for attenuation is commonly based on the standard‐linear‐solid (SLS) wave equation, in which case the quality factor (Q) has to be converted to stress and strain relaxation times. When using multiple attenuation mechanisms in the SLS method, it is difficult to directly estimate these relaxation time parameters. Based on a time domain complex‐valued viscoacoustic wave equation, we present an FWI framework for simultaneously estimating subsurface P wave velocity and attenuation distributions. Because Q is explicitly incorporated into the viscoacoustic wave equation, we directly derive P wave velocity and Q sensitivity kernels using the adjoint‐state method and simultaneously estimate their subsurface distributions. By analyzing the Gauss‐Newton Hessian, we observe strong interparameter crosstalk, especially the leakage from velocity to Q. We approximate the Hessian inverse using a preconditioned L‐BFGS method in viscoacoustic FWI, which enables us to successfully reduce interparameter crosstalk and produce accurate velocity and attenuation models. Numerical examples demonstrate the feasibility and robustness of the proposed method for simultaneously mapping complex velocity and Q distributions in the subsurface.

Multi-parameter wavefield reconstruction inversion for wavespeed and attenuation with bound constraints and total variation regularization

Wavefield reconstruction inversion (WRI) extends the search space of Full Waveform Inversion (FWI) by allowing for wave equation errors during wavefield reconstruction to match the data from the first iteration. Then, the wavespeeds are updated from the wavefields by minimizing the source residuals. Performing these two tasks in alternating mode breaks down the nonlinear FWI as a sequence of two linear subproblems, relaying on the bilinearity of the wave equation. We solve this biconvex optimization with the alternating-direction method of multipliers (ADMM) to cancel out efficiently the data and source residuals in iterations and stabilize the parameter estimation with appropriate regularizations.Here, we extend WRI to viscoacoustic media for attenuation imaging. Attenuation reconstruction is challenging because of the small imprint of attenuation in the data and the crosstalks with velocities. To address these issues, we recast the multivariate viscoacoustic WRI as a triconvex optimization and update wavefields, squared slowness, and attenuation factor in alternating mode at each WRI iteration. This requires to linearize the attenuation-estimation subproblem via an approximated trilinear viscoacoustic wave equation. The iterative defect correction embedded in ADMM corrects the errors generated by this linearization, while the operator splitting allows us to tailor ℓ1 regularization to each parameter class. A toy numerical example shows that these strategies mitigate crosstalk artifacts and noise from the attenuation reconstruction. A more realistic example representative of the North sea further shows the ability of the method to jointly reconstruct the wavespeed and the attenuation starting from a very crude attenuation-free initial model, the moderate strength of the crosstalk artifacts and the sensitivity of the multi-parameter reconstruction to noise.

Viscoacoustic least-squares migration with a blockwise Hessian matrix: an effective Q approach

Abstract We present an explicit inverse approach using a Hessian matrix for least-squares migration (LSM) with Q compensation. The scheme is developed by incorporating an effective Q-based solution of the viscoacoustic wave equation into a blockwise approximation to the Hessian in LSM, which is implemented after the so-called deabsorption prestack time migration (PSTM). The effective Q model used fully accounts for frequency-dependent traveltime and amplitude at the same imaging location. We can extract the effective Q parameters by scanning during previous deabsorption PSTM. This avoids the challenging task of building the Q model. The blockwise Hessian matrix approach decomposes the full Hessian matrix into a series of computationally tractable small-sized matrices using a localised approach. We derive the explicit formula of the offset-dependent Hessian matrix using an analytical Green's function obtained from deabsorption PSTM. In this way, we can approximate a reflectivity imaging for the targeted zone by a spatial deconvolution of the migrated result with an explicit inverse. The resulting scheme broadens the frequency-band of imaging by deabsorption, and improves the subsurface illumination and spatial resolution through the inverse Hessian. A high-resolution, true-amplitude migrated gather can then be obtained. Synthetic and field data sets demonstrate the proposed blockwise LS-QPSTM.

Fractional Laplacians viscoacoustic wavefield modeling with k-space-based time-stepping error compensating scheme

The spatial derivatives in decoupled fractional Laplacian (DFL) viscoacoustic and viscoelastic wave equations are the mixed-domain Laplacian operators. Using the approximation of the mixed-domain operators, the spatial derivatives can be calculated by using the Fourier pseudospectral (PS) method with barely spatial numerical dispersions, whereas the time derivative is often computed with the finite-difference (FD) method in second-order accuracy (referred to as the FD-PS scheme). The time-stepping errors caused by the FD discretization inevitably introduce the accumulative temporal dispersion during the wavefield extrapolation, especially for a long-time simulation. To eliminate the time-stepping errors, here, we adopted the [Formula: see text]-space concept in the numerical discretization of the DFL viscoacoustic wave equation. Different from existing [Formula: see text]-space methods, our [Formula: see text]-space method for DFL viscoacoustic wave equation contains two correction terms, which were designed to compensate for the time-stepping errors in the dispersion-dominated operator and loss-dominated operator, respectively. Using theoretical analyses and numerical experiments, we determine that our [Formula: see text]-space approach is superior to the traditional FD-PS scheme mainly in three aspects. First, our approach can effectively compensate for the time-stepping errors. Second, the stability condition is more relaxed, which makes the selection of sampling intervals more flexible. Finally, the [Formula: see text]-space approach allows us to conduct high-accuracy wavefield extrapolation with larger time steps. These features make our scheme suitable for seismic modeling and imaging problems.

An approach for attenuation-compensating multidimensional constant-Q viscoacoustic reverse time migration

Simulation of wave propagation in a constant-[Formula: see text] viscoacoustic medium is an important problem, for instance, within [Formula: see text]-compensated reverse time migration (RTM). Processes of attenuation and dispersion influence all aspects of seismic wave propagation, degrading the resolution of migrated images. To improve the image resolution, we have developed a new approach for the numerical solution of the viscoacoustic wave equation in the time domain and we developed an associated viscoacoustic RTM ([Formula: see text]-RTM) method. The main feature of the [Formula: see text]-RTM approach is compensation of attenuation effects in seismic images during migration by separation of amplitude attenuation and phase dispersion terms. Because of this separation, we are able to compensate the amplitude loss effect in isolation, the phase dispersion effect in isolation, or both effects concurrently. In the [Formula: see text]-RTM implementation, an attenuation-compensated operator is constructed by reversing the sign of the amplitude attenuation and a regularized viscoacoustic wave equation is invoked to eliminate high-frequency instabilities. The scheme is tested on a layered model and a modified acoustic Marmousi velocity model. We validate and examine the response of this approach by using it within an RTM scheme adjusted to compensate for attenuation. The amplitude loss in the wavefield at the source and receivers due to attenuation can be recovered by applying compensation operators on the measured receiver wavefield. Our 2D and 3D numerical tests focus on the amplitude recovery and resolution of the [Formula: see text]-RTM images as well as the interface locations. Improvements in all three of these features beneath highly attenuative layers are evident.

Viscoacoustic anisotropic wave equations

The wave equation plays a central role in seismic modeling, processing, imaging and inversion. Incorporating attenuation anisotropy into the acoustic anisotropic wave equations provides a choice for acoustic forward and inverse modeling in attenuating anisotropic media. However, the existing viscoacoustic anisotropic wave equations are obtained for a specified viscoacoustic model. We have developed a relatively general representation of the scalar and vector viscoacoustic wave equations for orthorhombic anisotropy. We also obtain the viscoacoustic wave equations for transverse isotropy as a special case. The viscoacoustic orthorhombic wave equations are flexible for multiple viscoacoustic models. We take into account the classic visocoacoustic models such as the Kelvin-Voigt, Maxwell, standard-linear-solid and Kjartansson models, and we derive the corresponding viscoacoustic wave equations in differential form. To analyze the wave propagation in viscoacoustic models, we derive the asymptotic point-source solution of the scalar wave equation. Numerical examples indicate a comparison of the acoustic waveforms excited by a point source in the viscoacoustic orthorhombic models and the corresponding nonattenuating model, and the effect of the attenuation anisotropy on the acoustic waveforms.

Modeling Frequency‐Independent Q Viscoacoustic Wave Propagation in Heterogeneous Media

AbstractQuantifying the attenuation of seismic waves propagating in the Earth interior is critical to study the subsurface structure. Previous studies have proposed fractional anelastic wave equations to model the frequency‐independent Q seismic wave propagation. Such wave equations involve fractional derivatives that pose computational challenges for the numerical schemes in terms of accuracy and efficiency when dealing with heterogeneous Earth media. To tackle these challenges, here we derive a new viscoacoustic wave equation, where the power terms of the fractional Laplacian operators are spatially independent, thus accurate and efficient methods (e.g., the Fourier pseudospectral method) can be adopted. Our derivation enables the resultant equation to capture both amplitude and phase signatures of the anelastic wave propagation by matching the complex wave numbers for all the frequencies of interest. We verify the derivation by comparing the dispersion curves of both the attenuation factor and the phase velocity produced by the new wave equation with their theoretical values as well as the Pierre Shale in situ measurements. Following that, we use a synthetic attenuating gas chimney model to demonstrate the attenuation effects on seismic waveforms and then construct the Q‐compensated reverse time migration to undo these effects for seismic image enhancement. Finally, we find that our forward modeling results can characterize the spatiotemporal attenuation effects revealed in the Frio‐II CO2 injection time‐lapse seismic monitoring data. We expect this proposed equation to be useful to quantify the attenuation in seismic data to push the resolution limits of seismic imaging and inversion.

Solving fractional Laplacian visco-acoustic wave equations on complex-geometry domains using Grünwald-formula based radial basis collocation method

We propose a radial basis function collocation method (RBF method) to solve fractional Laplacian visco-acoustic wave equation for the Earth media having heterogeneous velocity model and complex geometry. Unlike the fractional Laplacian wave equation proposed in Zhu and Harris (2014), the wave equation we consider has a different definition for the fractional Laplacian. Specifically, spectral and Riesz fractional Laplacians are considered in Zhu and Harris (2014) and the present paper, respectively. Accordingly, the Fourier pseudospectral method (FPS method) and the RBF method are employed to solve the spectral and the Riesz fractional Laplacian wave equations. The two wave equations are observed to produce obviously different wavefields. We demonstrate the validity and flexibility of the proposed RBF method by considering five benchmarks of seismic forward modeling: (1) two-dimensional Earth media with four types of velocity models (homogeneous, two-layer, homogeneous but complex-geometry, and heterogeneous models) and (2) a spherical medium with homogeneous velocity model. We make a three-way comparison among numerical solutions to the Riesz fractional Laplacian, the spectral fractional Laplacian, and the integer-order visco-acoustic wave equations, and observe that when wave attenuation is weak the Riesz wave equation yields more similar wavefield to that of the integer-order wave equation than the spectral wave equation does. Furthermore, uniform and quasi-uniform layouts for collocation points of the RBF method are considered, and the latter layout turns out to be economical since it can preserve the solution accuracy with the minimum number of collocation points. The RBF method is truly mesh-free and dimension-free and can easily handle high-dimensional, irregular domains. Additionally, the method is easier to implement than element-based methods, such as finite element and spectral element methods, for discretizing the Riesz fractional Laplacian.

Modeling the viscoelastic effects in P-waves with modified viscoacoustic wave propagation

Accurate full-waveform seismic modeling is a powerful tool for understanding wave propagation and building subsurface images. However, it can be computationally expensive for viscoelastic media. Viscoacoustic seismic modeling is much cheaper, but at the trade-off of using incomplete physics. We have developed a modified viscoacoustic wave simulation algorithm for modeling the viscoelastic effects of P-waves. The algorithm contains two viscoacoustic forward-modeling steps; the first is the same as the traditional viscoacoustic modeling, whereas the second propagation is generated using a residual error source, which is derived by comparing the viscoacoustic and viscoelastic wave equations in the form of stress-particle velocity formulations. The corrected P-wave particle velocities can be obtained by adding the wavefield from the second simulation step to the original (the first simulation step) viscoacoustic wavefield. Only P-waves are modeled. The overall cost is about twice that of viscoacoustic modeling, but it is significantly less than a viscoelastic propagation because there are fewer calculations, and we can use a coarser grid and larger time steps for the same accuracy. Numerical examples indicate that the P-wave waveforms, after correction, match those from viscoelastic wave modeling better than those from the original viscoacoustic simulation. Our method provides a cost-efficient alternative for approximating the viscoelastic effects in P-wave modeling.