Addressing viscous effects in acoustic full-waveform inversion

In conventional full-waveform inversion (FWI), viscous effects are typically neglected, and this is likely to adversely affect the recovery of P-wave velocity. We have developed a strategy to mitigate viscous effects based on the use of matching filters with the aim of improving the performance of acoustic FWI. The approach requires an approximate estimate of the intrinsic attenuation model, and it is one to three times more expensive than conventional acoustic FWI. First, we perform 2D synthetic tests to study the impact of viscoacoustic effects on the recorded wavefield and analyze how that affects the recovered velocity models after acoustic FWI. Then, we apply the current method on the generated data and determine that it mitigates viscous effects successfully even in the presence of noise. We find that having an approximate estimate for intrinsic attenuation, even when these effects are strong, leads to improvements in resolution and a more accurate recovery of the P-wave velocity. Then, we implement and develop our method on a 2D field data set using Gabor transforms to obtain an approximate intrinsic attenuation model and inversion frequencies of up to 24 Hz. The analysis of the results indicates that there is an improvement in terms of resolution and continuity of the layers on the recovered P-wave velocity model, leading to an improved flattening of gathers and a closer match of the inverted velocity model with the migrated seismic data.


INTRODUCTION
As seismic waves travel through the subsurface, the media in which they travel absorb part of their energy due to the inherent viscosity. These viscous effects are translated into an amplitude reduction and phase dispersion of the seismic waveforms (Causse et al., 1999).
Recorded data can be used in conventional full-waveform inversion (FWI) to recover high-resolution models of the physical properties by iteratively minimizing the L 2 -norm of the misfit between the recorded and the modeled data (Tarantola, 1984). To obtain quantitative models of the subsurface, viscous effects need to be considered in the inversion (Causse et al., 1999;Warner et al., 2012), but they are typically neglected in conventional 3D time-domain applications mainly due to the large computational cost of solving the viscoelastic wave equation and the necessity to use several relaxation mechanisms to have constant absorption with frequency (Blanch et al., 1995).
Seismic migration images are improved when the estimated velocity model is accurate and the viscous effects are accounted for, for example, in Q-compensated migration (Dai and West, 1994;Mittet et al., 1995;Causse and Ursin, 2000;Zhu et al., 2014), which requires an estimate for intrinsic attenuation, described by the quality factor Q. In the past few decades, FWI has had some success in obtaining high-resolution and accurate velocity models, as well as smooth intrinsic attenuation models of the subsurface, despite its ill-posed nature. With this aim, several authors have implemented viscoacoustic FWI in the frequency domain, in which viscous effects are incorporated in seismic modeling by introducing intrinsic attenuation as the imaginary component of P-wave velocity (Pratt, 1990(Pratt, , 1999Pratt and Worthington, 1990;Song et al., 1995;Plessix et al., 2012;Prieux et al., 2013;Operto et al., 2015), whereas others have implemented FWI in the time domain (Cheng et al., 2015;Plessix et al., 2016;Bai et al., 2017). A variety of approaches have been used to recover the intrinsic attenuation, including compensating the gradient for viscous effects (Xue et al., 2016(Xue et al., , 2017, running tomography for Q after having performed tomography and FWI for velocity only (Zhou et al., 2013), and using a semiglobal inversion algorithm . Kurzmann et al. (2013) recently apply a passive viscoacoustic FWI strategy in which intrinsic attenuation was incorporated into the modeling as a fixed parameter in a time-domain inversion. This resulted in better P-wave velocity models of the subsurface, but required modeling of viscous effects in the time domain at each iteration of the inversion, which increased the computation cost when compared with conventional acoustic FWI.
Here, we implement an alternative scheme by adapting the method of Agudo et al. (2016Agudo et al. ( , 2017, originally designed to mitigate elastic effects, to address viscous effects. The aim of this modified method is not to obtain an intrinsic attenuation model of the subsurface, but to mitigate viscous effects and improve the estimated P-wave velocity models from FWI. Unlike the original approach for elastic effects, the method introduced herein requires modeling the viscous effects in the time domain for all shots before carrying out the inversion. Our approach also requires an approximate estimate of Q for P-waves (Q p ), here referred to as Q for simplicity, which is efficiently obtained for field data using Gabor transforms of near-offset data (Wang, 2004) prior (or simultaneously) to the application of the workflow. As a result, the current method provides a data set in which the viscous effects are mitigated. Analogously to the method for elastic effects, the present method also involves the computation of matching filters. Nevertheless, their characteristics are different to those computed for elastic effects (Agudo et al., 2017) given the different characteristics of each type of medium and, consequently, different characteristics in wave propagation. The inversion of the resulting data is then performed using the acoustic wave equation without including viscous effects, and we show that this results in better quality of the recovered models than those obtained after acoustic FWI of the observed data. Additionally, the present approach reduces the computational cost of the forward modeling during the inversions. Although the present approach has the potential to address the viscous and elastic effects at the same time, here we focus on mitigating viscoacoustic effects and we demonstrate its application on the synthetic and field data.
In terms of the structure of this paper, we first show the impact of not considering the viscous effects in the recovered P-wave velocity models after acoustic FWI of 2D synthetic models. Then, we demonstrate how the new approach for incorporating viscous effects effectively improves the estimation of the velocity models from data affected by intrinsic attenuation with synthetic and real-data exam-ples. Finally, we discuss the limitations and further applications of this method.

METHODOLOGY
The amplitude and phase of the recorded seismic data are affected by intrinsic attenuation (Causse et al., 1999). This is depicted in Figure 1, which compares the recorded pressure signal at 1490 m from a source in a homogeneous medium with constant velocity (V P ¼ 2000 m∕s) and with or without viscosity. P-waves are shown for an acoustic (the black curve) and a viscoacoustic medium using a quality factor of Q ¼ 20 (the gray curve). Due to viscous effects, the signal generated in a viscoacoustic medium has lower amplitudes than the signal in an acoustic medium and is delayed in time (note the difference in time of the peaks and troughs of both signals). These effects are often ignored in acoustic FWI, which can lead to inaccurate estimates of P-wave velocity models.
To account for these effects and obtain improved recovered P-wave velocity models of the subsurface, herein we introduce a method for mitigating viscous effects in FWI using matching filters. The proposed method consists of suppressing viscoacoustic effects from the observed data by computing matching filters that match modeled viscoacoustic to modeled acoustic data generated from estimates of the subsurface P-wave velocity and Q models, followed by the application of the filters to the observed data. This is analogous to the method proposed by Agudo et al. (2016Agudo et al. ( , 2017 to mitigate elastic effects in acoustic FWI, but with three main differences: (1) It requires estimation of Q, which can be performed simultaneously and independently to an acoustic inversion of the observed data and does not need to be updated at each iteration of the workflow, (2) it does not require an estimate of V P ∕V S that depends on a recovered V P model, and (3) viscoacoustic data instead of elastic data are modeled, in which the characteristics of wave propagation are different. Given that the viscoacoustic wave propagation is dispersive, the spectrum changes with time and stationary filters cannot capture this variation. Thus, it is necessary to account for the nonstationary behavior of viscoacoustic wave propagation, whereas this is not the case when mitigating elastic effects.
The first difference in the methodology with respect to that of Agudo et al. (2016Agudo et al. ( , 2017 is a consequence of having to model viscoacoustic wave propagation. For this purpose, an estimated Q model is required, which can be obtained independently from an inversion. Several methods have been reported in the literature to estimate Q from reflection seismic data (e.g., Quan and Harris, 1997;Zhang and Ulrych, 2002;Wang, 2004). Here, we use the method of Wang (2004) to estimate Q for surface field data because it provides reliable and robust trace-by-trace estimations when compared with other methods (Lupinacci and Oliveira, 2015).
The second difference is due to the fact that elastic effects are not considered and also because an approximate Q model can be obtained directly from reflection data without having to perform any inversion. We demonstrate that such an approximate Q model is enough to mitigate most of the viscoacoustic effects from the data, leading to improved P-wave velocity models.
Finally, the nonstationary behavior in the recorded data due to intrinsic attenuation needs to be accounted for, especially when the viscoacoustic effects are strong given the stationary nature of the matching filters. This is done by introducing a modified workflow that is more robust against strong viscoacoustic effects and also by using short overlapping filters. Figure 1. Pressure signal recorded at a certain distance from the source in a homogeneous medium with a constant P-wave velocity of 2000 m∕s when there is no intrinsic attenuation (black line) and when there is intrinsic attenuation (Q p ¼ 20, gray line), which delays and attenuates the recorded signal.

Data-matching workflows
To further explain these differences, we now introduce the workflow in Figure 2 referred to as workflow 1and to explain the steps in more detail, we present two modified workflows to address strong viscoacoustic effects and strong noise in the dataworkflows 2 and 3, respectivelyand we explain how overlapping filters are computed.
The first step of workflow 1 (step 1) consists of performing conventional acoustic FWI, as introduced by Tarantola (1984) (any other functional can be used in principle, such as, for example, the amplitude-balanced functional of Shen, 2010). This preliminary estimate of the velocity model is affected by the intrinsic attenuation because it does not consider the viscous effects; however, it is more accurate than the starting model.
A model of Q is then simultaneously estimated from the data, and the estimated Q model and the recovered P-wave velocity model are then used to generate acoustic and viscoacoustic modeled data (step 2). The latter is more accurate than the former, i.e., it is closer to the observed data, because it accounts for the dispersion and attenuation due to the viscoacoustic effects.
Then, we use these modeled data for correcting the observed data affected by intrinsic attenuation. This correction is carried out with matching filters. We compute the matching filters that map the viscoacoustic into the acoustic data and then we convolve them with the observed data, leading to a corrected (or matched) data set (step 3). This data set contains the same events as the observed data, but with a difference in terms of amplitude and phase related to the suppression of viscoacoustic effects, which depends on the accuracy of the estimates of Q and V P .
At this stage, we perform acoustic FWI of the matched data set using the original starting model, which leads to an improved P-wave velocity model given that the inverted data set is less affected by viscoacoustic effects than the original observed data (step 4).
To further improve the quality of the recovered velocity model, the workflow can be repeated again by generating acoustic and viscoacoustic data with the new recovered V P model and the estimate of Q. This process can be repeated several times until a target level of accuracy is achieved (step 5)for instance, the workflow can be repeated until a certain root-mean-square (rms) difference of the matched data with respect to the observed data for all shot gathers is achieved, although this is not further investigated here. Nevertheless, tests on synthetic and field data lead to the conclusion that the first iteration of the workflow has the largest effect on mitigating viscous effects because it is the first time that intrinsic attenuation is not ignored.
In the following sections, we also validate two variations of workflow 1, namely, workflows 2 and 3, which improve the performance of our method in the presence of strong absorption effects and random noise in the observed data, respectively, without any added costand even a lower cost for workflow 2.
Workflow 2 consists of eliminating the first acoustic inversion in step 1 of the workflow in Figure 2 and performing the remaining steps of the workflow by using the starting model instead of the estimated velocity model as well as a Q estimate obtained from the data. In particular, we estimate Q by using Gabor transforms of reflection data (Wang, 2004). Viscoacoustic effects are then mitigated by using the estimated Q model and the starting P-wave velocity model directly in steps 2−5 of the workflow in Figure 2; i.e., acoustic and viscoacoustic data are modeled for all the shots prior to any inversion, matching filters that match the viscoacoustic to the acoustic data are computed, these are applied to the observed data, and the resulting matched data set is inverted using acoustic FWI. Hence, the cost of workflow 2 is slightly higher than that of a conventional acoustic inversion and, if another iteration of the workflow in step 5 is implemented, the computation time is similar to that of one single datamatching workflow iteration because the number of acoustic inversions is also two, but the viscoacoustic effects are mitigated twice.
Workflow 2 is suitable in the presence of strong viscous effects, in which the recovered P-wave velocity model obtained after step 1 of the workflow may contain artifacts due to crosstalk between the P-wave velocity and intrinsic attenuation. Instead, in workflow 2, acoustic FWI is performed only after having mitigated the viscoacoustic effects. In fact, a Q model is estimated directly from the data, and acoustic and viscoacoustic data are modeled with different propagators prior to any inversion, thus reducing any possible correlation between the P-wave velocity and Q models. Thus, an acoustic inversion of the resulting matched data results in a recovered velocity model that is less affected by crosstalk effects because most viscoacoustic effects have been mitigated prior to it.
On the other hand, step 1 is performed normally in workflow 3, and in step 4, a starting model for the acoustic inversion is obtained Figure 2. Proposed workflow used to mitigate viscous effects (referred to in the text as workflow 1).
by heavily smoothing the model recovered after the first acoustic inversion in step 1 of the workflow. This helps in mitigating the impact of noise on the recovered velocity models because the starting model used in step 4 is in principle more accurate than the original starting modelprovided the viscoacoustic effects are not very strongand also because any rapid variations in the recovered velocity model obtained after step 1 due to noise in the data are smoothed.
We validate the performance of workflows 1-3 on the synthetic examples, and we apply workflow 2 on a field data set because the level of noise is low compared with the signal. The extension of the current method to address viscoelastic effects is straightforward and would require estimates of Q p , Q s , and V P ∕V S to model viscoelastic data in step 2 of the workflow.
The cost of workflow 1 introduced herein is two to three times the cost of a conventional acoustic FWI, mainly due to the cost of the additional acoustic inversions in steps 1 and 4. This is, however, reduced when implementing workflow 2, so that the total elapsed time is slightly above the compute time of a single conventional acoustic inversion. As a result, workflow 2 has a lower cost than a viscoacoustic FWI, given that the latter requires viscoacoustic wave modeling for each iteration in the inversionin the time domain, this requires considering several relaxation mechanisms, thus increasing the computation time. The present method is therefore an intermediate solution between acoustic and viscoacoustic FWI because it can partially account for viscous effects without inverting for Q at a cost that is lower than a full viscoacoustic inversion in the time domain.

Data-matching strategy and parameters
As in the method described by Agudo et al. (2016Agudo et al. ( , 2017, matching filters are treated as optimum Wiener (1949) filters. We estimate these filters from the least-squares minimization of the following quadratic error function: where w is a single Wiener filter that matches M traces simultaneously (M is an odd number and is our first matching parameter), N is the number of time samples per trace, and P is the order of the filter. The quadratic function in equation 1 matches the modeled viscoacoustic data (d va ) to the modeled acoustic data (d ac ) and imposes smoothness of the matched data to avoid artifacts caused by trace-by-trace matching. As in the original derivation (Agudo et al., 2016(Agudo et al., , 2017, this equation assumes small variations between neighboring traces and introduces smoothing on the matched data set. Minimization of the error function with respect to the filter coefficients yields a set of normal equations that involve crosscorrelations and autocorrelations of both modeled data, which are solved using Levinson's (1947) recursion method given the symmetry of the equations.
Our strategy is to compute matching filters on a trace-by-trace basis close to the edges of the data set and increase the number of simultaneously matched traces linearly as we move away from the first or last traces of each shot gather. Instead of using single filters with the length of the trace, we use overlapping short filters weighted using a Blackman window (Chakraborty, 2013) such that the sample at the center of the filter has the largest weight. This reduces the computation costs significantly and mitigates artifacts caused by matching entire traces that have different spectra by smoothing the matched data along each trace.
The use of short Wiener filters is also adequate because it partially accounts for the nonstationary behavior of the recorded data due to intrinsic attenuation because each filter in equation 1 is stationary. Thus, a combination of overlapping short filters in time ensures that the matched data are only locally stationary. The length of the filters is, hence, our second matching parameter. The computation of the matched data is carried out in parallel, distributing the tasks across multiple nodes. The parallel estimation of the Wiener filters results in a compute time of approximately less than 10% (depending on the length of the filters) the computing time of conventional acoustic FWI for the same data set.
We want to highlight that this approach is different to that of Warner and Guasch (2016) and that of Zhu and Fomel (2016), in which matching filters that match the modeled to the observed data are estimated during the inversion, and these are adapted as the model changes for each iteration to derive the earth parameters. In these approaches, the matching filter converges to a Dirac delta function when the estimated model converges to the reference model. Instead, here we estimate the impact of viscoacoustic effects on the data, and mitigate the latter using matching filters prior to any inversion. Thus, the matching filters are not used here to derive subsurface parameters, but to suppress the effect of a physical property from the data.
In the following section, we test our approach on synthetic and field data and we discuss the optimal choice of the data-matching parameters.

RESULTS
We now demonstrate the benefits of the present method on synthetic and field data. First, we discuss the impact of different parameters, e.g., the length of the filters. The use of trace normalization on the FWI functional during the inversion is also investigated on a simple layered 2D model. Viscoacoustic effects on a more realistic 2D synthetic data set are then examined and mitigated using matching filters. We also investigate the impact of the presence of noise in the data when using the method outlined herein. In the final example, we apply the suggested strategy to a field data set and analyze the significant improvements obtained after addressing viscoacoustic effects on this data set, which may lead to easier and more accurate interpretations of the subsurface structure. For all tests in this manuscript, density is related to velocity following Gardner et al.'s (1974) law.

2D horizontally layered models
We apply our method to horizontally layered 2D models of the subsurface. Figure 3 shows the true P-wave velocity model, which consists of two constant-velocity layers with a change in velocity at the sea bottom (600 m depth) and at 850 m depth, and which is used to generate acoustic data. The same figure also shows two different intrinsic attenuation models that are used to generate viscoacoustic data: Model A consists of two layers with constant and relatively low-quality factors, whereas the quality factor in model B changes gradually with depth and it has larger values. Thus, the data modeled with the intrinsic attenuation model A will contain stronger viscous effects (as absorption is inversely proportional to the quality factor), and will be used to assess the impact of strong intrinsic attenuation on acoustic FWI and also on the performance of the current method. On the other hand, the intrinsic attenuation model B simulates a different environment with smaller viscous effects. For both models, we use a large value for the quality factor within the water layer to simulate a marine data set, given absorption is negligible in water.
First, we generate acoustic and viscoacoustic data using the 2D laterally invariant models in Figure 3. These models are 15 km long by 1.2 km deep. The data are generated with different acoustic and viscoacoustic modeling implementations that have the same order of accuracy (second order in time and sixth order in space). The grid spacing is 10 m, which ensures that there is little dispersion of P-waves. The time dependency of the source is defined with a Ricker wavelet with peak frequency at 15 Hz. The records are generated with a time sampling of 1 ms and 3 s of time length. We generate data setting off 41 sources at 6 m depth uniformly distributed across the model. The records of the pressure signal are sampled with 967 receivers at 16 m depth, which are also uniformly distributed across the model. To emulate a marine-recording setting, we impose a free-surface boundary condition at the top of the model. This guarantees that surface-related multiples, as well as source and receiver ghosts, are present in the generated data. At the remaining lateral and bottom boundaries, we impose absorbing layers, preventing any spurious reflections to be reinjected into the simulation domain.
The generated acoustic and viscoacoustic data are then inverted with conventional acoustic FWI using a multiscale approach so that the data are inverted from 3 up to 15 Hz in a total of 60 iterations and using the starting velocity model in Figure 3, which leads to the recovered velocity models shown in Figure 4a obtained at the center of the 2D models. We observe that the model after acoustic FWI of acoustic data matches the true model closely: This is the reference result because viscous effects are neglected in the modeling and inversion. On the other hand, neglecting intrinsic attenuation during acoustic FWI of the viscoacoustic data leads to recovered velocity models that are slower on average and less accurate, especially with increasing depth. This is an effect of amplitude attenuation and time dispersion introduced by absorption of seismic waves in the subsurface. As expected, the impact of viscous effects is larger for model A (the red line in Figure 4a) due to the smaller values for the quality factor of the reference model when compared with those in model B.
Next, we use the recovered models after acoustic FWI of viscoacoustic data to mitigate viscoacoustic effects using the suggested workflow 1 in Figure 2. For each of the recovered models, we generate acoustic and viscoacoustic data separately and we use a rough estimate of the Q models by averaging the true intrinsic attenuation models within the subsurface, as depicted by the dashed red and blue lines in Figure 3.
Then, we find the matching filters that match the viscoacoustic to the acoustic modeled data, and we apply these filters to the corresponding observed (viscoacoustic) data. We design the matching filters in this section to be 325 ms in length, and we match 13 traces simultaneously so that the spacing between the first and the last simultaneously matched trace is 180 m. The dashed red line and solid blue line in Figure 4b show the results of performing acoustic FWI on the respective matched data sets for models A and B. We observe that using the current workflow and average estimates of Q is sufficient to significantly mitigate viscoacoustic effects on the recovered velocity models for both data sets. The recovered models are now faster and more accurate than those in Figure 4a before the workflow was applied; i.e., they are closer to the reference model. However, the result for model A (the dashed red line) still shows a residual imprint caused by the larger intrinsic attenuation effects.
To further reduce viscoacoustic effects for model A, we perform a second iteration of the workflow in Figure 2 by using the newly recovered velocity model and the same average Q estimate to generate acoustic and viscoacoustic data. This leads to the recovered model shown in Figure 4b (the solid red line) after using the original starting model for the acoustic inversion. The recovered model in Figure 4b shows that a second iteration of the workflow produces a recovered model that is closer to the reference model when viscoacoustic effects are large, but it requires a longer computation time as a second workflow iteration means a second acoustic inversion.
To minimize the computation time and still obtain a good recovered model in the presence of strong viscoacoustic effects, we now apply workflow 2, which consists of eliminating the inversion in step 1 of the workflow in Figure 2 and using the starting model to mitigate viscoacoustic effects, initially. If we apply workflow 2 to the viscoacoustic data for model A and we use the original starting model and the average Q model to mitigate viscoacoustic effects, we obtain the recovered velocity model represented with the solid cyan line in Figure 4b. This model is more accurate than that obtained with workflow 1 after a single iteration (the dashed red line) at approximately the same cost, and it is very close to the model obtained after two iterations of workflow 1 (the solid red line), even though the cost of workflow 2 is much lower because it requires one less acoustic inversion. Thus, we conclude that this is the best approach for mitigating large viscous effects when using the method outlined in this paper. Figure 3. Vertical profiles of the reference P-wave velocity and Pwave quality factor (solid lines), the starting velocity model (dashed black line) and the averages of the Q models (color dashed lines) for the two horizontally layered models.
Alternatively to the workflow introduced herein, Shen (2010) proposes the use of an amplitude-balanced FWI functional to mitigate viscoacoustic (and elastic) effects. The aim of this functional is to enhance phase rather than amplitude fitting of the observed and modeled data. Figure 4c and 4d shows the inversion results when carrying out acoustic FWI with the amplitude-balanced functional and using our workflow 1. Figure 4c depicts the result of acoustic FWI of the observed acoustic and viscoacoustic data, whereas Figure 4d shows the result of acoustic FWI of the matched data after one iteration of workflow 1 for intrinsic attenuation models A and B note that we do not need to apply workflow 2 because the amplitude-balanced approach reduces some of the viscous effects.
The results in Figure 4c are more accurate than those in Figure 4a, except for the acoustic data because we are disregarding amplitude information generated with the same modeling code as for the inversion, but there is still an imprint of viscous effects in depth and at the sea bottom, especially for model A. Thus, we apply the workflow in Figure 2 (i.e., workflow 1) to mitigate intrinsic attenuation effects for both viscoacoustic data sets by using the recovered velocity models in Figure 4c and the average Q models in Figure 3. Thus, we generate acoustic and viscoacoustic data, derive matching filters, compute the matched data, and perform acoustic FWI of the matched data. This results in the recovered models in Figure 4d after just one iteration of the workflow. We observe that the resulting velocity models are more accurate than those in Figure 4c before addressing viscous effects. This improvement is also observed at the sea bottom. Then, we conclude that the amplitude-balanced functional of Shen (2010) has significant benefits when dealing with the mitigation of viscoacoustic effects in FWI, but it might not be sufficient in the presence of strong viscoacoustic effects, in which phase differences introduced by intrinsic attenuation are important. This is because the latter does not accurately account for viscoacoustic wave propagation. Instead, we perform acoustic and viscoacoustic modeling in step 2 of the workflow introduced herein, thus accounting for larger amplitude and phase variations. Moreover, the proposed data-matching workflow can be combined with any functional, e.g., the amplitude-balanced functional of Shen (2010), to further mitigate the viscoacoustic effects, as shown in the previous example. Hence, in the remaining tests in this manuscript, we use an amplitude-balanced functional combined with any of our three proposed workflows.
The accuracy of the recovered velocity models using our proposed method also marginally depends on the choice of the matching parameters, i.e., the length of the filters and the number of traces matched simultaneously. We now assess the quality of the matched data with respect to the choice of the matching parameters for the viscoacoustic data set for model B by comparing the matched data with the observed acoustic and viscoacoustic data. Figure 5 compares the observed viscoacoustic data and the matched data for different matching parameters with the acoustic observed data. We only show the effect of the matching parameters on model B here, but similar conclusions can be drawn for the same tests on model A. This is expected because the velocity model does not change, and only the intrinsic attenuation is different. In fact, differences with respect to the acoustic data are slightly larger for model A, as the viscoacoustic effects are stronger, but these do not alter the conclusions obtained for the tests performed on model B.
In Figure 5a, we display the viscoacoustic data for model B, whereas the difference of the latter with the acoustic data free of viscous effects is presented in Figure 5b. We observe that the difference is small for the direct arrivals that have traveled through the water layer. This is expected because the quality factor Q is very high in this layer; hence, the wavefield does not get attenuated. The same is observed for refracted arrivals. This difference is very significant for reflections and deeper refractions due to the time delay introduced by viscous effects with respect to the events in the acoustic data set free of viscous effects and also due to amplitude changes.
On the other hand, data differences with respect to the acoustic data are reduced after the application of the current method. Figure 5c shows the matched data obtained after step 3 of the workflow in Figure 2, in which we match data trace-by-trace and we use overlapping short filters with a length of 325 ms, whereas its difference with the true acoustic data is displayed in Figure 5d. By comparing Figure 5b and 5d, we observe a clear overall decrease in the magnitude of the data difference after the application of the current approach, which is mainly due to a phase correction and is in part due to an amplitude adjustment; e.g., note the amplitude correction of the refraction at the second interface (the black arrow in Figure 5b when compared with the same event in Figure 5a). Nevertheless, we observe that there is a small distortion of the events in which different waves interfere, such as the area inside the black circle in Figure 5c. In these regions, datamatching trace-by-trace leads to spurious effects because the spectrum of the acoustic and viscoacoustic traces is very different, and for this reason they are mismatched.
To reduce these effects and impose continuity and smoothness of the data, we match multiple traces simultaneously; i.e., we look for the Wiener filter that optimally matches viscoacoustic to acoustic data for a selection of neighboring traces and a given time window. Figure 5e shows the matched data obtained after simultaneously matching 13 traces so that the spacing between the first and the last matched trace is of 180 m, whereas the length of the filters is unchanged. This results in a smoother matched data set, in which data-matching artifacts are reduced: We note there are no spurious effects due to data matching inside the circle in Figure 5e. After several experimental tests on multiple data sets, we concluded that the simultaneous matching of traces in a range between 150 and 250 m leads to the best compromise between not having artifacts in the estimated filters and generating a filter that does not predict the viscoacoustic effects.
Finally, we study the impact of the length of the overlapping filters on the resulting matched data. To test its impact, we compute matched data by matching 13 traces simultaneously and with lengths of the Wiener filters set to 50 and 1500 ms. The resulting matched data with 50 ms filters' length are shown in Figure 5g, and its difference with respect to the true acoustic data is presented in Figure 5h. Figure 5i displays the matched data obtained with filters of 1500 ms length, whereas Figure 5j depicts the difference of the latter with respect to the true acoustic data. Comparing Figure 5g and 5h and Figure 5i and 5j with those in Figure 5e and 5f for an intermediate length of the filters, we observe 1) Shorter filters lead to slightly noisier data and larger differences for reflections (black arrow on the right side in Figure 5h) and refractions (black arrow on the left side in Figure 5h).  Figure 3 for the (a) observed viscoacoustic data and the matched data obtained (c) trace by trace and with a length of the filters of 325 ms, and the matched data obtained by matching 13 traces simultaneously and with filter lengths of (e) 325, (g) 50, and (i) 1500 ms. The right column shows the difference of the data on the left column with respect to the true acoustic data, generated without intrinsic attenuation.
2) Longer filters lead to smoother and more continuous data, but with larger errors close to the source (the horizontal black arrow in Figure 5j) and also at events at late times (the black arrows on the left side in Figure 5j).
The difference in the performance of the different filter parameters displayed in Figure 5 can be better understood by comparing a single trace in each of the previous data sets with that in the acoustic data set. Thus, Figure 6 compares each trace for the data sets in Figure 5 (the gray lines) with the same trace in the acoustic data (the black line) at 2 km positionat a relatively long offset. Figure 6a shows the expected reduction in amplitude and phase delay introduced by viscoacoustic effects (the gray line) when compared with the acoustic trace (black line). The differences in amplitude and time delay are mitigated when applying the proposed strategy, and this can be observed in Figure 6b-6e. Figure 6b displays the extracted trace for the matched data set obtained by matching trace by trace and with a length of the filters of 325 ms (as in Figure 5c). There is no time delay in the matched trace in Figure 6b, which will help the inversion to converge to a global minimum, and amplitude differences are now minimal. Figure 6c shows the comparison for the matched data set (gray line) obtained with the same length of the filters but matching 13 traces simultaneously (as in Figure 5e). Although the resulting matched trace is very similar to that in Figure 6b, matching multiple traces simultaneously improves the performance of the filters in regions where multiple events cross each other (see the circled regions in Figure 5c and 5e). Figure 6d and 6e compares the acoustic data (the black line) with the matched data (the gray lines) obtained by matching 13 traces simultaneously and shorter (50 ms) or longer (1500 ms) lengths of the filters to those used in Figure 6c, respectively. From these, we observe that shorter filters correct for the phase difference but perform less well in matching amplitude differences at early times, whereas longer filters lead to an improved amplitude correction for first arrivals, but these deteriorate for late times.
Because FWI mainly uses wide-angle refracted arrivals to build a velocity model (Pratt et al., 1996), we conclude that long-matching filters (longer than half a second) are better suited than short filters (less than 0.1 s) and that filters of intermediate length (0.2-0.5 s) produce an optimum match for early and late arrivals.
The latter is due to the fact that long matching filters do not account for the nonstationary behavior of recorded data due to intrinsic attenuation, whereas intermediate filters partially account for it because they are only applied to a reduced part of the data in time.
To further illustrate this, Figure 7 compares several matching filters that lead to the matched traces in Figure 6c and 6e, obtained for lengths of 325 and 1500 ms. More specifically, Figure 7a shows the acoustic and viscoacoustic data obtained after step 2 of workflow 1 for model set B corresponding to the matched data in Additionally, the choice of filtering parameters also has an impact on the time required to compute the matched data. This is illustrated in Table 1, which shows a comparison between the computation times required to obtain the previous matched data sets for all the sources and receivers in the data set using a single node with 24 cores. The latter shows that the length of the filter has the largest contribution to the elapsed timewe recall that the compute time for the Levinson's recursion method is proportional to the square of the number of samples in the data being filteredand that increasing lateral smoothing also increases the computation time.
Thus, after this test and the other tests in this manuscript, we recommend the use of filters with an intermediate length between 300 and 400 ms, which provide a good compromise between computation time and quality of the matched data.  Figure 5) between the true acoustic data (black) and the viscoacoustic observed data (gray) and (b-e) between the true acoustic data (black) and the matched data (gray) for the shot gathers in Figure 5c-5j. The matched data in (b-e) correspond to that obtained (b) trace by trace and with length of the filters of 325 ms and that obtained by matching 13 traces simultaneously with lengths of the filters of (c) 325, (d) 50, and (e) 1500 ms. The matched trace in (b and c) is closest to the true acoustic data.

2D synthetic tests
We now validate our approach on a more realistic subsurface data set obtained from the marine Marmousi2 model (Martin et al., 2006), and we assess the effect of noise present in the observed data on the final recovered velocity model. We also compare the improved recovered models with those obtained with viscoacoustic FWI of the observed viscoacoustic data keeping the model of intrinsic attenuation fixed during the inversion.
We first generate observed acoustic and viscoacoustic data from the original P-wave velocity model ( Figure 8a) and a smooth intrinsic attenuation model, with values of the quality factor that range from 20 at the sea bottom to 100 at the bottom of the model (Figure 8c). The water layer is 450 m deep and we generate an acoustic signal for 99 sources uniformly spaced every 160 m, and placed at 6 m depth. The sources are set off sequentially. The observed data for each source are recorded by 791 receivers uniformly distributed every 20 m. Acoustic and viscoacoustic data are generated with different numerical implementations of the constitutive laws. However, they have the same order of accuracy in space and time. The source wavelet is given by a Ricker wavelet with a peak at 11 Hz. The useful frequencies for inversion are in the range of 2−20 Hz. The record length is 6 s, and the time-sampling interval is 1 ms. The true model is discretized with 10 m grid spacing along the horizontal and vertical directions. The discretization parameters guarantee that there is no grid dispersion of acoustic waves. We apply absorbing boundaries on the sides and on the bottom of the model. At the top of the model, we apply a free-surface boundary condition, which adds surfacerelated multiples as well as source and receiver ghosts, mimicking these events when carrying out marine seismic surveys.
Next, we perform acoustic FWI of the acoustic and viscoacoustic observed data by using the amplitude balanced functional (Shen, 2010) because this functional is more robust when dealing with viscous effects, as pointed out in the previous section. We use a multiscale approach by inverting the data from 3 up to 18 Hz and use the smooth starting model depicted in Figure 8b obtained by strongly smoothing the true P-wave velocity model in Figure 8a. Figure 8e and 8f shows the recovered velocity models after acoustic FWI of the acoustic and viscoacoustic observed data, respectively. The high-fidelity and accurate velocity model in Figure 8e is the reference solution. In this case, the constitutive law used for generating the data matches the one that is used for carrying out the inversion. Hence, the main factors limiting the accuracy of the inverted velocity model are the band-limited nature of the signal and the acquisition geometry.
On the other hand, the recovered model in Figure 8f illustrates that having an inaccurate constitutive law for carrying out the inversion, in this case neglecting viscoacoustic effects, has an adverse impact on the recovered P-wave velocity model. This result also demonstrates that having a more robust functional (amplitudebalanced) does not suffice to compensate inaccurate physics in the inversion. In attenuating media, the energy envelope propagates more slowly and the seismic events are recorded at later times. Table 1. Elapsed compute times to obtain the matched data sets in Figure 5c, 5e, 5g, and 5i. The distance between receivers is 15 m (i.e., 13 traces represent 180 m distance from first to last matched traces). Long filters and large lateral smoothing increase the computation time.  Figure 6c and 6e, respectively. Short overlapping filters account for the nonstationary behavior of intrinsic attenuation.
In addition, because the physics is incomplete, the waveforms have the wrong phase due to differences in the constructive and destructive interference. In summary, ignoring viscous effects when estimating the wavefield leads to an inverted velocity model that is slower on average; it is less sharp and lacks some of the structures observed in the model in Figure 8e. We further assess the quality of the recovered models by computing the relative error for each point within the dotted box in these figures with respect to the true model, and we display the average value on the top right corner, which gives an indication of the overall model misfit within this region. We note that this value is lower in Figure 8e with respect to that in the starting model in Figure 8b, and it is in between these two values for the recovered model in Figure 8f. Consequently, we apply the proposed workflow in Figure 2 (i.e., workflow 1) to further mitigate viscoacoustic effects. We use an estimate of the Q model that is a two-layer average of the true Q model and which is more accurate at the top, as displayed in Figure 8d, and we model acoustic and viscoacoustic data using this  estimate and the recovered model after acoustic FWI of viscoacoustic data (Figure 8f). We then compute matching filters of 350 ms length that match the latter to the former modeled data by matching 11 traces simultaneously so that the spacing between the first and the last matched trace is 200 m. These filters are then applied to the observed (viscoacoustic) data, and the resulting matched data are inverted again with acoustic FWI and the same original starting model, which results in the recovered velocity model in Figure 8g.
Despite the use of a very approximate Q model, the data correction introduced by workflow 1 with a single iteration of the proposed workflow results in a velocity model that is faster, more accuratenote the smaller average relative errorand with a higher level of detail than that in Figure 8f. As a comparison, we perform viscoacoustic FWI of the observed viscoacoustic data by using the same average Q model during the inversion and keeping the latter fixed, and the recovered velocity model is shown in Figure 8h. This results in a similar velocity model to that in Figure 8g. This model, which is used here as a reference model, is slightly better resolved given that we have used the same constitutive laws for generating and inverting the data, whereas different constitutive lawsviscoacoustic and acoustic wave equations, respectivelywere used for modeling and inversion leading to the velocity model in Figure 8g. Nevertheless, we note that the data-matching method recovers a more accurate and slightly higher velocity of the horizontal layer at approximately 2.1 km depth on the left side of the model because this is thicker in the true model.
The differences for the velocity models in Figure 8 at two horizontal positions of 4 and 8 km are illustrated in Figure 9a and 9b, respectively. The latter shows the improvement in accuracy obtained by acoustically inverting the matched data (the blue line with circles) rather than the viscoacoustic data (the red line with crosses)for instance, at 2.3 km depth in Figure 9a, or at 1.6 km depth in Figure 9b. These well-log comparisons also show that the recovered model obtained after acoustic FWI of matched viscoacoustic data (the blue line with circles) is closer to the result of viscoacoustic FWI of viscoacoustic data (the orange line with squares) than that obtained without mitigating viscous effects (the red line with crosses).
In practice, there are three advantages of using the current method with respect to a full viscoacoustic inversion with a fixed Q model in the time domain: 1) All inversions with the current method are acoustic, which reduces the computation timei.e., we do not need to model several relaxation mechanisms at each iteration of the inversion to have constant intrinsic attenuation with frequency when modeling in the time domain. 2) It can be combined with the method described by Agudo et al. (2016Agudo et al. ( , 2017 to account for viscoelastic effects. 3) It provides a matched data set, in which viscous effects are mitigated for all arrivals.  Finally, Figure 10 shows a comparison between the observed viscoacoustic data for a shot at the center of the model (Figure 10a) and the matched viscoacoustic data used in the previous tests (Figure 10b). Their difference with respect to the observed acoustic data is shown in Figure 10c and 10d, respectively. Comparing these figures, we observe that there is a visible correction in terms of amplitude for all the events; the correction in phase is less obvious from these plots because the data set is now more complicated than the 2D horizontally layered models. As expected, the overall data differences with respect to the observed acoustic data, affected by intrinsic attenuation, are reduced in Figure 10d with respect to Figure 10c due to the correction for viscoacoustic effects, especially in the middle region of the data set.
This reduction can also be quantified in terms of the average rms difference of the data in Figure 10a and 10b with respect to the acoustic data: A value of 0.16 is obtained for the viscoacoustic data in Figure 10a and 10c, whereas a value of 0.12 is obtained for the matched viscoacoustic data in Figure 10b and 10d, which is approximately a 25% difference. Acoustic FWI of the matched data then leads to the recovered model in Figure 8g. Thus, as a result of applying workflow 1 introduced herein on the observed data, one can observe a significant improvement in the inverted model in Figure 8g.

Impact of noise
We now analyze the effect of adding strong random noise on the observed viscoacoustic data and repeating the previous tests with the noisy viscoacoustic data. As a result, we apply workflow 3 to assess its performance in mitigating the impact of noise on the recovered models. Figure 11a shows a shot in the center of the model for the viscoacoustic observed data after adding strong random noise. The latter is added by computing the mean of the absolute difference between minimum and maximum values for each trace in a shot gather and adding a random percentage of the latter that is between 0% and 50% to every sample in every shot gather. Additionally, the noise is filtered with the spectrum of the source wavelet to suppress noise with frequencies outside the bandwidth of the data. Figure 11b shows the matched observed data after applying the workflow in Figure 2 (workflow 1) once and using the same inversion and modeling parameters as in the previous tests, but now the velocity model used for the modeling is different. As in Figure 10, we observe an amplitude correction of all the events, but the amplitude of noise is also modified, which may negatively affect the quality of the recovered models. Additionally, we want to highlight that we do not apply a mute on the matched data note that there are no data before a certain time for each tracebut we compute matching filters only within the region in which the modeled acoustic data in step 2 of the workflow are nonzero in order to speed-up the computation. Figure 11c and 11d displays the difference of the panels above with the same gather for the noise-free acoustic data. These data differences also show that there is a correction in terms of amplitude and phase for the reflections and refractions, but the noise increases this difference in some regions of the data. To determine whether the corrections in amplitude and phase of the main events have more importance than the boosting of noise, we perform acoustic FWI of the noisy viscoacoustic data and acoustic FWI of the matched noisy viscoacoustic data after one iteration of the workflow using the original starting velocity model in Figure 8b, and the results are shown in Figure 12a and 12b. Strong noise clearly introduces artifacts into the recovered models when compared with the noise-free results (Figure 8f and 8g), which is particularly detrimental on the result of acoustic FWI of the matched data (Figure 12b). Despite the overall improvement in terms of velocity accuracy in some regionswe note an overall increase in velocity, for instance, on the horizontal layer at approximately 2 km depth and from 2 to 8 km distancestrong noise introduces artifacts in certain areas of the model, and it damages the shape of the dome at the bottom of the model. Thus, we apply workflow 3, which is useful when the level of noise is high. Figure 11. Representative shot gathers at the center of the Marmousi2 model for (a) the viscoacoustic data with noise and (b) the matched viscoacoustic data obtained from the noisy observed data by simultaneously matching 11 traces and with a length of the filters of 350 ms. Their difference with respect to the true acoustic data is shown in the bottom row.
We perform the acoustic inversion in step 4 of the workflow with a starting model that is built by heavily smoothing the recovered model in Figure 12a (after step 1 of the workflow), which results in the starting model in Figure 12c. Running our workflow with this starting model leads to an improved recovered velocity model as depicted in Figure 12d. This inverted velocity model is faster on average and more accurate. Note the decrease in the average relative error with respect to Figure 12b and the higher sharpness of the model when compared with the one obtained if viscoacoustic effects are neglected during the inversion (Figure 12a). This result is comparable with that obtained after applying workflow 1 to the noise-free observed data (Figure 8g).
On the other hand, the test in Figure 12e shows that performing acoustic FWI on the observed viscoacoustic data, starting from the model in Figure 12c does not lead to an inversion result as good as the one in Figure 12d, where the matched data were used. Hence, workflow 3 is well suited and robust when the observed data are strongly contaminated with noise.
As a final comparison with a reference solution, Figure 12f shows the recovered model obtained from viscoacoustic FWI of the noisy viscoacoustic data using the original starting model and keeping Q fixed during the inversion. The model of intrinsic attenuation used is equal to the estimated Q model in Figure 8d, resulting in a comparable velocity model. We note again that this has been obtained using the same constitutive laws for generating the observed data and inverting it, and that time-domain viscoacoustic FWI takes longer than time-domain acoustic FWI because the constitutive law includes several relaxation mechanisms, which need to be included when carrying out wavefield simulations.

Field data
The proposed methodology to mitigate viscous effects is now implemented on a 2D seismic line acquired in the Carnarvon Basin, on the Australian North West Shelf (Kalinicheva et al., 2017). This survey was designed and optimized with the view of inverting the Figure 12. Vertical slices of the recovered P-wave velocity models obtained after acoustic FWI of (a) the noisy viscoacoustic data and (b) the matched viscoacoustic data inverted with the original starting model, (c) a new starting model obtained by smoothing the model in (a) and the recovered velocity models obtained from acoustic FWI using the starting model in (c) for (d) the matched data and (e) the noisy viscoacoustic data. (f) The recovered model after viscoacoustic FWI of the noisy viscoacoustic data using the original starting model. acquired data with FWI. In this survey, a single low-frequency source with a peak frequency at 8 Hz and useful signal at less than 2.5 Hz was towed at 10 m depth to inject a pressure signal into the water, and the shot spacing was 50 m.
The data were recorded with a 10,000 m long streamer cable towed at 25 m depth with a receiver spacing of 12.5 m. The acquisition was done in a marine environment in which water depth ranged from 600 to 1600 m. Thus, these data contain strong water-bottom multiples and significant refracted energy due to the long-recorded offsets, and they also have a good signal-to-noise ratio at low frequencies.
Our implementation of VTI anisotropic acoustic FWIwith fixed anisotropic parameters, not shown here for simplicityinvolves minimal data preprocessing: The data are band-pass filtered between 1.5 and 24.5 Hz, and a source wavelet is extracted from the data by using matching filters and the direct arrivals present in the data. To mitigate viscoacoustic effects, we implement the current method. First, we estimate an intrinsic attenuation model of the subsurface by using the method of Wang (2004), which is based on Gabor transforms of reflection data. In this method, the Gabor transform of the data is first computed, which reveals the time-varying frequency characteristics of each trace. Assuming plane-wave propagation, each transformed reflection trace is then fitted to an intrinsic attenuation in the least-squares sense, which results in an estimate of Q.
To apply it to a field data set that is not a pure reflection data set, we extract the nearest-offset trace for each shot and we use this to estimate Qthe separation between each shot and the first receiver in the streamer is approximately 100 m. We then correct the reflection data for spherical divergence using an rms velocity model built from the contractor's starting velocity model and we use the resulting data, shown in Figure 13a, as an input to estimate Q. The latter is estimated on a trace-by-trace basis in the time domain for nine different time layers, and the obtained Q model in time is converted to depth using the starting velocity model and a standard time-to-depth conversion workflowthe latter uses the starting model of interval velocities in depth, it converts the latter to average velocities in time, and it transforms the Q model in time to depth.
Finally, the Q model in time is extended laterally to cover the entire length of the velocity model and Gaussian smoothing across three cells is applied eight times to avoid sharp discontinuities, which results in the Q model in Figure 13b. This is a smooth and approximate intrinsic attenuation model with values generally increasing with depth. Due to the characteristics of this setting (e.g., no gas cloud or fine layering), we assume viscous effects are larger at shallow depths because the sedimentary rocks are less compacted and, hence, more energy is lost due to absorption. The values of Q range from 20 near the sea bottom up to 420 at 4 km depth.  We want to stress that this model is approximate for the following reasons: (1) The data used to estimate Q using Gabor transforms are not a pure reflection data set, which is required in the method proposed by Wang (2004) that has been used here, but a data set containing the nearest offset traces, (2) surface-related multiples are kept in the data because removing them could adversely modify the spectra of the traces and thus the estimation of Q, (3) an approximate rms velocity model is used to correct for spherical divergence and convert the Q model from time to depth, and (4) the Q model is estimated trace by trace using a limited number of time layers to reduce the computation time. Nevertheless, such an approximate Q model should be sufficient to mitigate intrinsic attenuation given the results from our synthetic examples.
The strategy used to mitigate viscoacoustic effects using matching filters for this data set is based on the synthetic data results, the estimated Q model, and the good signal-to-noise ratio of the recorded data: We use a single iteration of workflow 2 to account for possible strong intrinsic attenuation. We do not use workflow 3 given the low level of noise present in the data, and, hence, we use the original starting model in all acoustic inversions. Thus, the computational cost is slightly over the cost of a single acoustic inversion, and less than the cost of two acoustic inversions. Modeling of acoustic and viscoacoustic data is performed without considering anisotropythe model of Q is isotropicbut all inversions for this data set are anisotropic, given that velocity models are anisotropic (the models of anisotropic parameters used are those provided by the contractor). Figure 14a and 14b depicts the recorded data for two different shots, showing a very good signal-to-noise ratio. These shots were acquired along the same sailing line with a distance separation of 22.5 km. We note the level of noise is very low before the first arrivals and also later in timeit is only noticeable within the first traces of the second shot, close to the source, which is coherent noise. The application of the current strategy to mitigate viscoacoustic effects results in the two matched shot gathers in Figure 14c and 14d, in which we have matched 15 traces simultaneously so that the spacing between the first and the last matched trace is 175 m, and the length of the filters is set to 325 ms. We observe there are few visible artifacts due to data matching given our choice of matching parameters. By comparing Figure 14c and 14d with Figure 14a and 14b, respectively, we observe that there is a visible correction of the amplitude; for instance, compare the amplitude change of the refracted arrivals inside the rectangular box, and some of the noise in Figure 14b for traces close to the source is mitigated in Figure 14d by the smoothing introduced by matching multiple traces simultaneously. It is, however, more difficult to observe phase changes at this scale because the data set contains multiple events and changes are of the order of milliseconds.
The recovered models after anisotropic acoustic FWI of the observed and matched data sets reveal the benefits of using the current strategy to mitigate viscoacoustic effects. The inversions are carried out in a multiscale fashion iterating in blocks of four iterations for each frequency band. The first frequency band is limited at 2.5 Hz with a cut-off filter. The frequency band is then progressively widened, after the completion of each block of iterations, up to 24 Hz. The starting velocity model is based on prestack time migration stacking velocities (Kalinicheva et al., 2017) and is displayed in Figure 15a. Performing acoustic FWI of the observed data results in the velocity model in Figure 15b, in which different layers and features not present in the starting model can now be identified; for instance, we note the presence of channels in the shallow part of the model, a major unconformity that crosses the model between 2 and 2.5 km depthdenoted by a black arrow in Figure 15b and also strong and irregular reflectors at the bottom of the model with high velocities at less than 2.5 km depth, which are interpreted to be Figure 15. Vertical slices for the field data set of (a) the starting velocity model, (b) the recovered velocity model after anisotropic acoustic FWI of the observed data, and (c) the recovered velocity model after anisotropic acoustic FWI of the matched observed data.  Figure 15b and 15c, respectively. The color scale is the same as in Figure 15.
Addressing viscous effects in FWI R625 basaltic intrusions. Finally, the recovered model after acoustic FWI of the matched data using workflow 2i.e., we mitigate viscoacoustic effects directly from the starting modelis displayed in Figure 15c. When compared with the model in Figure 15b, mitigating viscoacoustic effects with the current strategy leads to a model with higher velocities on average, especially below the unconformity, in which the layers are better resolvede.g., note the thinning of the fast layer right above the unconformityand in which structures below the unconformity can be better identified and followed across the model.
We further compare the two recovered velocity models by magnifying within the area inside the green box in Figure 15c, and these are respectively displayed in Figure 16a and 16b. In the section in Figure 16b, we clearly observe the presence of high-velocity layers and clear vertical discontinuities within these layers (the blue arrows in Figure 16b)the latter might indicate the presence of faultsbelow the unconformity between 10 and 20 km distance, which are difficult to observe in Figure 16a before accounting for viscous effects.
To assess the quality of the recovered model in Figure 15c (or Figure 16b) and analyze whether it represents an improvement with respect to that obtained after conventional acoustic FWI, we perform 2D anisotropic Kirchhoff poststack depth migration. First, we migrate the observed data using the recovered model after acoustic FWI of the observed data in Figure 15b and we show the corresponding surface-offset common-image gathers (CIGs) in Figure 17a. Then, we migrate the Q-compensated observed datausing the Q model in Figure 13b and band-pass filtering the resulting data with the same filter that was applied on the observed datawith the velocity model obtained after acoustic FWI of the matched data, and the corresponding CIGs are displayed in Figure 17b. We observe that the CIGs are slightly flatter in Figure 17b with respect to Figure 17a, and there is amplitude equalization with depth and thinning of the events after accounting for viscous effects. We note that a more sophisticated migration algorithm that considered intrinsic attenuation during migration could be applied instead of applying Q compensation, but we do not have access to such software.
Finally, Figure 18a and 18b shows the migrated stacked sections respectively obtained from the recovered velocity models after acoustic FWI of the observed data and after acoustic FWI of the matched data, shown in black and white, and the corresponding velocity models in Figure 15b and 15c are overlaid (in color). The section in Figure 18b in comparison with that in Figure 18a shows a higher level of detail and better continuity of the layers, especially below the unconformity. Furthermore, the layers we observed in the velocity model below the unconformity in Figure 16b are also present in the migrated section in Figure 18b between common depth points 4000 and 5000and they are in the same location. We conclude that the proposed strategy is effective for mitigating viscous effects in FWI when inverting data affected Figure 17. Common offset gathers for the observed data (a) before and (b) after Q compensation with the Q model in Figure 13b and after band-pass filtering, obtained by using the velocity models in Figure 15b and 15c, respectively. Figure 18. Stacked migrated sections (gray color scale) obtained from the observed data (a) before and (b) after Q compensation and using the recovered velocity models obtained from acoustic FWI of the observed and the matched data, respectively. The velocity models in Figure 15b and 15c are overlaid. by intrinsic attenuation, leading to an improvement of the inverted velocity models.

DISCUSSION
We have shown that the suggested data-matching workflow can be used to address viscoacoustic effects in (anisotropic) acoustic FWI on synthetic and field data sets. Unlike inverse Q-filtering, our approach can address intrinsic attenuation effects and remove them from the data without the 1D-layered earth assumption (Wang and Guo, 2004) and, hence, it can handle real geologic complexity. The results presented in this manuscript demonstrate the benefits of this approach for velocity model building. We have also introduced and validated a modified workflow to account for strong intrinsic attenuation effects and noisy data (even for low levels of coherent noise, as was the case for the field data set). This requires an initial matching step after generating data from the starting model directly and using a starting model obtained by smoothing the recovered model after acoustic FWI of the observed data. The tests performed on the synthetic data sets also demonstrated that an accurate estimate of a smooth background Q model is enough to mitigate most of the viscoacoustic effects. This has been also confirmed on a field data set in which an estimate of the intrinsic attenuation model has been obtained from near-offset data without removing multiples from the data. The latter has been avoided because it may modify the spectra of the data and degrade the estimate of Q.
Additionally, we also suggest that the current workflow can be extended to include viscoelastic effects in acoustic FWI by combining the current method and that suggested by Agudo et al. (2016Agudo et al. ( , 2017 to account for elastic effects. A generalized data-matching workflow based on that in Figure 2 would involve the modeling of acoustic and viscoelastic data in step 2 and would require an estimate of Q for P-and S-waves as well as an estimate of V P ∕V S . With this aim, we propose the use of a local/global hybrid method to estimate Q p , Q s , and V P ∕V S from the observed data in a sparse grid similar to the approach presented by da . In this method, velocity and intrinsic attenuation are simultaneously inverted using a combination of local gradient descent methods to update velocity and a global approach to generate random and smooth models of intrinsic attenuationmore specifically, da Silva et al. (2017) suggest using a quantum particle swarm optimization algorithm. The latter could potentially be used to obtain smooth models of the above parameters needed to mitigate viscoelastic effects with the current method, thus yielding a matched data set with mitigated viscoelastic effects, followed by conventional acoustic FWI. We expect that this should result in an improved P-wave velocity model at a cost of between two to three times the cost of a single acoustic inversion, well below the cost of viscoelastic inversion because the latter requires elastic modeling, which increases the compute time significantly (Chapman et al., 2010).
Future research will therefore focus on the application of this methodology to address viscoacoustic and viscoelastic effects in three dimensions.

CONCLUSION
Due to absorption in the subsurface, seismic waves are attenuated and dispersed and this is recorded at the surface. If these phenomena are not included in the constitutive law when carrying out FWI, the wavefield is then computed with errors and these translate into inaccuracies in the inverted velocity models. In this paper, we have extended a methodology based on matching filters, which was originally developed to mitigate elastic effects, to address viscous effects. We have demonstrated that our workflow leads to improvements in recovered models of P-wave velocity after acoustic FWI. We have first analyzed the impact of viscoacoustic effects on the recovered P-wave velocity models, and we have then shown a substantial improvement in terms of resolution and accuracy. This improvement is observed especially in the deeper parts of the inverted models, after the application of the current method on the horizontally layered and laterally variant 2D synthetic data sets. These results also indicate that an approximate estimate of Q is enough to mitigate most of the viscoacoustic effects. We have also suggested two modified workflows that perform well when there is strong intrinsic attenuation and when the observed data are contaminated with noise. These alternative workflows have a computational cost from one to slightly over two times that of acoustic FWI. The application to a field data set demonstrates that our approach has potential on a realistic case, in which Q is estimated from near-offset data using Gabor transforms. When compared with the result of anisotropic acoustic FWI, layers deep in the model that are difficult to track in the former are better imaged using the current workflow. This is verified by improvements in the flattening of CIGs and a migrated section that is easier to interpret due to improved resolution and continuity of events. Application of the current method to 3D data sets is straightforward: It requires data modeling in 3D and data matching by carefully considering the acquisition geometry. Future tests will focus on applying this methodology to account for viscoelastic effects in three dimensions.