Massively parallel frequency-domain FWI (FFWI code)

We implemented a 3D massively parallel frequency-domain FWI code (called FFWI ) with the goal to perform the first application of frequency-domain FWI on a real industrial OBC dataset when the forward problem is solved with sparse direct solver . Frequency-domain FWI is more specifically suitable for stationary-recording acquisitions such as OBC or OBN because they provide a broad angular illumination of the subsurface which allows us to limit FWI to a limited number of frequencies (see the resolution analysis of FWI in the tutorial page). In the frequency domain, wavefield simulation is a boundary value problem which requires to solve a large and sparse linear system with multiple right-hand sides. Two categories of linear algebra methods can be used: direct or iterative methods.The advantage of direct methods is to give the solution in a finite number of operations (this is particularly beneficial for ill-conditioned Helmholtz problems for which iterative solvers need to be equipped with efficient preconditioner) and the efficient processing of multiple right-hand sides (RHSs) by forward/backward substitution. The drawback is to involve a RHS-independent "preprocessing" step (the LU decomposition) which is memory demanding and computationally expensive. However, the time complexity of one LU decomposition \((O(N^6))\) (namely, for one frequency) is the same as that of an explicit time stepping method for a dense surface seismic acquisition, where \(N\) stands for the number of points along one dimension of a \(N^3\) grid and the number of sources scales to \(O(N^2)\). Moreover, the memory demand of one LU decomposition \((O(N^4))\) is one order of magnitude smaller than the memory storage of a time-domain dataset \((O(N^5))\). Accordingly, there was no obvious theoretical reason to prevent the assessment of this technology for 3D stationary-recording acquisition such as OBC or OBN, in particular considering recent advances carried out by the community developing sparse direct solver to decrease the complexity of the problem with low-rank compression strategies (e.g., Amestoy et al., 2016).

Figures 1-4 illustrate some results of 3D visco-acoustic VTI frequency-domain FWI on Ocean Bottom Cable (OBC) data collected in the North Sea. Both seismic modelling and inversion have been performed in the frequency domain with FFWI . The linear system resulting from the discretization of the time-harmonic wave equation has been solved with the finite-difference method of Operto et al., 2014 and the sparse multifrontal direct solver MUMPS using quite limited computational resources provided by the computer center SIGAMM hosted by Observatory of Côte d'Azur. This frugal use of computational resources may sound counter intuitive as sparse direct solvers are considered as being computionally expensive to tackle 3D problems due to the memory demand of LU factorization and its limited scalability for large scale problems. In contrast, we have shown with the MUMPS team the efficiency of this approach to tackle problems involving a few tens millions of unknowns, hence validating our feasibility analysis published in 2007 ( Operto et al., 2007). See Operto et al., 2014; Operto et al., 2015; Amestoy et al., 2016; Operto & Miniussi, 2018 for more details.

Figure 1: Seismic imaging of an oil field in the North Sea by frequency-domain FWI. (Left): P-wave velocity model of the subsurface across a low-velocity gas cloud. (Right): A 10Hz monochromatic wavefield computed by solving the Helmholtz equation with the sparse multifrontal direct solver MUMPS is superimposed. A total of 685 computer cores are typically used to perform such simulations and FWI at 10Hz. ( Operto et al., 2015; Amestoy et al., 2016; Operto & Miniussi, 2018).

Figure 2: On the resolution power of FWI. The above figure shows some horizontal and vertical sections of P-wave velocity models of the oil field. (Left) Velocity model obtained by reflection traveltime tomography. This model was used as starting model for FWI. (Right) Velocity model obtained by frequency-domain FWI. (a-c) the slices cross-cut (a) sand channel deposits at 175m depth, (b) scarves left by drifting icebergs on the paleo sea bed at 500m depth, and (c) a gas cloud at 1km depth. The two bottom panels show vertical sections across (d) the gas cloud and (e) its periphery ( Operto et al., 2015; Amestoy et al., 2016; Operto & Miniussi, 2018).

Figure 3: (a) Velocity model of the oil reservoir. (b) Q model of the oil reservoir reconstructed by frequency-domain FWI. The imprint of attenuation in the seismic data fit is highlighted in Figure 4. Overall, the values of Q are consistent with the expected geology with a positive correlation between low velocities in the soft sediments and the gas cloud and high attenuating zones ( Operto & Miniussi, 2018).

Figure 4: Comparison between recorded and simulated data. The synthetic seismograms have been simulated with the same frequency-domain modeling engine as that used to perform FWI before inverse Fourier transform. The recorded data are plotted with a red/white/blue color scale. The simulated data are superimposed with black variable area. The two sets of synthetic are in phase if the black overprint the red. The simulated data are computed in FWI models when (a) attenuation has not be taken into account during FWI and (b) attenuation is updated during FWI. We show that dispersive waves generated by a shallow wave guide are better matched when attenuation is taken into account (yellow arrow) ( Operto & Miniussi, 2018).