^{1}

^{2}

^{1}

^{*}

^{1}

^{2}

Edited by: Catherine McGeoch, D-Wave Systems, Canada

Reviewed by: Henri Calandra, Total, France; Amandeep Singh Bhatia, Purdue University, United States; Matthias Möller, Delft University of Technology, Netherlands; Travis S. Humble, Oak Ridge National Laboratory (DOE), United States

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Quantum computing exploits quantum mechanics to perform certain computations more efficiently than classical computers. Current quantum computers have performed carefully tailored computational tasks that would be difficult or impossible for even the fastest supercomputers in the world. This “quantum supremacy” result demonstrates that quantum computing is more powerful than classical computing in some computational regimes. At present, it is unknown if any computational problems related to the Earth's subsurface fall within these regimes. Here, we describe an approach to performing seismic inverse analysis that combines a type of quantum computer called a quantum annealer with classical computing. This approach improves upon past work on applying quantum computing to the subsurface (via subsurface hydrology) in two ways. First, the seismic inverse problem enables better performance from the quantum annealer because of the Earth's relatively narrow distribution of P-wave velocities compared to the broad distribution of hydraulic conductivities. Second, we develop an iterative approach to quantum-computational inverse analysis, which works with a realistic set of observations. By contrast, the previous method used an inverse method that depended on an impractically dense set of observations. In combination, these two advances significantly narrow the gap a quantum-computational advantage for a practical subsurface geoscience problem. Closing the gap completely requires more work, but has the potential to dramatically accelerate inverse analyses for subsurface geoscience.

Computation has played a critical role in subsurface geoscience for decades. It has been used to simulate ocean circulation (Pinardi et al.,

Recent computer performance trends indicate that performance improvements are diminishing, suggesting that the rising tide of improving computational performance may be coming to an end. This trend has led to increased use of novel computational methods such as graphical processing unit computing (Fatemi and Poppe,

The improved scaling behaviors of certain quantum algorithms is what made possible the recent demonstration of quantum supremacy, where a quantum computer performed a calculation that would push the fastest classical supercomputers beyond their limits (Arute et al.,

This work, which does not attempt to show any quantum advantage, explores two applications—seismic inverse analysis and hydrologic inverse analysis. It improves upon previous work by enabling the solution of more realistic problems. Our approach enables 2D seismic inverse analysis, whereas previous work focused on a 1D, layered approach (Souza et al.,

This work also goes a step beyond the empirical observation of the performance and identifies problem characteristics that enable better performance for the quantum computer. These insights can help guide future work to find an advantage for quantum computers in applications to the subsurface. However, since quantum annealing hardware is still in its early stages, the problems we look at are still relatively simple compared to similar problems solved using classical computing techniques. These problems are scaled according to current quantum-computing hardware's computational size, and advances in hardware will allow for more complex problems to be addressed.

If the long-term goals of this research are successful, subsurface geoscience will be able to exploit the theoretical advances that have been demonstrated for calibrating models to data (Wiebe et al.,

The remainder of this manuscript is organized as follows. In Section (2), we describe quantum annealing, which is the quantum algorithm that we leverage, and our approach to formulating inverse analysis as a problem suitable for quantum annealing. Section 3 describes the results of applying this approach to seismic and hydrologic problems. A discussion of various aspects of our approach, including a problem characteristic that improves performance and possibilities with future quantum hardware and methods is presented in Section 4. Finally, concluding remarks are made in Section 5.

Quantum annealing is a heuristic optimization algorithm, similar to simulated annealing (Kirkpatrick et al.,

The input to a D-Wave quantum annealer is a vector _{i}) and a matrix _{ij}). The matrix,

From a practical perspective, the quantum annealer can be thought of as minimizing a function of the form

where each spin, _{i}, is either −1 of +1, and is called the Ising formulation. On the D-Wave annealer, there are additional sparsity constraints on the matrix _{ij}. That is, some of the _{ij} are constrained to be zero. Further, the D-Wave hardware limits the range of values—called the dynamic range—that can be set for _{i} and _{ij}, and problems with coefficient outside these ranges have to be rescaled to fit. This rescaling can have the effect of pushing small coefficients into the hardware's noise range. A more accurate description of the behavior of the annealer is that it is drawing from a distribution that preferentially samples values of

where each bit, _{i}, is either 0 or 1, and is related to the formulation in Equation 1 via _{i} = 2_{i} − 1. This change, from inverting values of _{i} ∈ {−1, 1} as in Equation 1 to _{i} ∈ {0, 1} as in Equation 2, is what differentiates the Ising model from the QUBO model. The values of _{i} and _{ij} can also be transformed into the values of _{i} and _{ij} by associating like terms in Equations (1) and (2). Equations (1) and (2) are two equivalent ways of formulating the same problem. We will use the formulation in Equation 2 throughout.

We consider an inverse approach where the goal is to divide the subsurface into two different materials with constitutive properties based on measurements obtained on the boundary of the domain, which aligns well with the quantum annealer's ability to perform binary optimization. Using standard variable names from seismic inversion, the objective function used in the inverse analysis takes the form

where _{low} or _{high}, which will correspond to the low and high values of the relevant parameter field (P-wave velocity for the seismic problem, or hydraulic conductivity for the hydrology problem). We attempt to minimize this global objective function using an iterative process, where each iteration involves solving a related optimization problem that is suited to the quantum annealer. We begin with an initial guess, ^{(0)}, and iteratively produce ^{(k+1)} from ^{(k)}.

During an iteration, we find the model update by creating a quadratic unconstrained binary optimization (QUBO) problem that approximates Equation (3). This QUBO problem can then be solved with the quantum annealer. The QUBO objective function takes the form

where ^{(k)}. This estimate is given by

where _{m} is defined to be

where _{i} is the ^{th} standard basis vector of size ^{th} component. Essentially, _{m} is an operator that flips the value of its input and then forward models it. This makes the computational cost of each iteration equal to the cost of _{i} = 0 then _{i} = 1 then _{i} = 1 then

The least-squares objective function for iteration

Note that this equation has a constant term, a term that is linear in ^{M}, is used to update the model, ^{(k+1)}, using Equation (7). We used default values for the annealing time, thermalization time, and post-processing. D-Wave's heuristic embedder was used to embed the problem graph on the D-Wave, which generally involves embedding a complete graph.

We study the steady-state groundwater flow equation,

where _{h} or _{l} where _{h} is a high conductivity value and _{l} is a low conductivity value. The inverse analysis's goal is, given a set of hydraulic head observations, to infer the spatially variable hydraulic conductivity, _{h} or _{l}. This process corresponds to determining where two different materials exist in the aquifer. If the highly conductive material were sand and the low conductivity material were clay, the inverse analysis would answer the questions: where is the sand? and where is the clay?

To obtain a hydrologic inverse problem, we generate two hydraulic conductivity fields based on two real-world examples where the hydraulic conductivity at each location is either _{h} or _{l}. We first look at the example from Lu and Robinson (

We have a 700 × 700 meter domain in this application, with seven discrete permeability blocks in the _{l} or _{h}. Given this problem geometry, there are 2^{49} ≈ 5.6 × 10^{14} possible solutions. The exact model we use is shown in

The two hydrology permeability models used in this experiment, where yellow locations are high permeability and green locations are for low permeability values. We place seven discrete permeability block locations in the x-direction and seven in the y-direction. The domain is 700 × 700 m, so each discrete permeability block is 100 × 100 m. The computational mesh grid is dx = dy = 10 meters. We have 24 receivers, denoted by red triangles, which are spread in a checkerboard pattern across the top of the domain.

In this application, our goal is to find the distribution of P-wave velocity values that give rise to a set of wavefield measurements from receivers on the surface of the domain. We study the acoustic wave equation,

where _{high} or _{low}, so this problem can be thought of as locating two different materials in the subsurface. We look at two different examples: one is a salt body in an constant background medium, and the other is a two-layer faulted example. In this application, our domain of interest is 1 kilometer wide and 0.5 km deep and includes 50 possible velocity value locations: 5 in the vertical direction and 10 in the horizontal direction. To keep computational costs similar to that of the hydrologic inverse analysis, our experiment only uses one source, and we choose our source ^{50} ≈ 1.1 × 10^{15} possible solutions.

The seismic velocity models used in this experiment, where yellow locations are high velocity and green locations are for low velocity values. We place ten discrete velocity block locations in the

We applied the approach to inverse analysis previously discussed to seismic and hydrologic problems. For each of these physical problems, we consider a case where _{high}/_{low} is large and another where _{high}/_{low} is relatively small. For the seismic problem, the case where _{high}/_{low} is relatively small is realistic. On the other hand, the case where _{high}/_{low} is large is realistic for the hydrologic problem.

_{0} = _{low}. The stopping criteria for iterations is when the same model output was selected for two iterations in a row. There is a strong tendency to converge to a good result when the contrast is low. In the low contrast setting, the inverse approach gets all the hydraulic conductivities correct in four analyses and at most 5 incorrect in the remaining six analyses. On the other hand, there is a strong tendency to converge to a lackluster result in the high contrast case. In the high contrast setting, the inverse approach gets all the hydraulic conductivities correct three times, but gets 12 or more incorrect in the remaining seven analyses. The performance in these two settings indicates that the high contrast case is more challenging for the inverse method than the low contrast case.

The convergence of our quantum annealing inverse approach for 10 low _{0} = _{low}. The numbers at the end convergence point represent the number of model runs that converged to that model value. The approach works well in the low contrast regime, but there are problems in the high contrast regime due to the quantum annealer's limited dynamic range. Realistic hydrologic inverse problems have a high contrast, so this is problematic.

The convergence of our quantum annealing inverse approach for 10 low _{0} = _{low}. The numbers at the end convergence point represent the number of model runs that converged to that model value. The approach works well in the low contrast regime, but there are problems in the high contrast regime due to the quantum annealer's limited dynamic range. Realistic seismic inverse problems have a low contrast, so these problems are well-suited to the quantum annealer.

For the low contrast case, we use the same velocity values and contrast as used in the initial example: _{low} = 4, 250_{high} = 4, 750_{low} = 2, 000_{high} = 5, 000_{0} = _{low}. The stopping criteria for iterations is when the same model output was selected for two iterations in a row. Like the hydrologic inverse analysis, the low contrast problem shows better convergence behavior than the high contrast problem. In the low contrast setting, the inverse approach gets all the velocities correct five times, between one and four velocities incorrect four times, and seven velocities incorrect once. Again, the performance of the inverse approach declines as the contrast increases. In the high contrast setting, the inverse analysis never obtains all the velocities correctly. The number of incorrect velocity values tends to cluster around seven, which was the worst result in the low contrast case. Due to the inconsistent nature of the convergence pattern in low vs. high contrast cases and seismic vs. hydrologic examples, we choose to not provide an algebraic form for the convergence patterns.

Inverse analysis in subsurface flow problems is challenging for a variety of reasons. One source of challenges is the high variability in hydraulic conductivity and permeability found in the Earth's subsurface. ^{−7} milliDarcy to 10^{4} milliDarcy (Lie,

Variability in permeability in unconsolidated sediments (Fetter,

Clay | 10^{−11}–10^{−8} |

Silt, sandy silts, clayey sands, till | 10^{−8}–10^{−6} |

Silty sands, fine sands | 10^{−7}–10^{−5} |

Well-sorted sands, glacial outwash | 10^{−5}–10^{−3} |

Well-sorted gravel | 10^{−2}–10^{−4} |

Note the extreme variability between sediment types (nine orders of magnitude) and the variability within sediment types (typically two orders of magnitude).

This extreme variability adds to the challenges for the quantum annealer because high contrasts in the parameters result in more variability in the QUBO coefficients.

The cumulative distribution function for an Ising problem in

While different geologic units may have large contrasts in hydraulic conductivity, there is much less variability in P-wave velocity values both between and within rock types than in hydraulic conductivity, as seen in

Variability in P-wave velocity in common rock types in exploration seismology (Bourbié et al.,

Wet sands | 1, 500 – 2, 000 |

Saturated shales and clays | 1, 100 – 2, 500 |

Porous and saturated sandstones | 2, 000 – 3, 500 |

Limestones | 3, 500 – 6, 000 |

Salt | 4, 500 – 5, 500 |

Note that there is much less variability in P-wave velocity values both between and within rock types than in Hydraulic Conductivity of unconsolidated sediments in

One of the most notable limitations of the current work is that the resolution of the subsurface parameterization is limited. The hydrologic inverse problem had a 7 × 7 grid of hydraulic conductivities, and the seismic problem had a 5 × 10 grid of P-wave velocities. This work was done with D-Wave's older generation 2000Q hardware. D-Wave's current generation Advantage hardware uses a Pegasus graph (Dattani et al.,

It should also be noted that further methodological developments could improve the resolution of the inverse model. The full domain could be explored by moving through the domain in a tiling fashion. Another possibility would be to use an alternative to Equation (4) that has some natural sparsity. For example, parameters associated with regions that are physically distant from each other might tend to have a small quadratic term, which could be neglected in some cases. Many possibilities cannot be explored here – this is just the beginning of using quantum computing for subsurface applications.

We have considered the application of noisy, intermediate-scale quantum computing to subsurface geoscience. In particular, we have used a quantum annealer to solve seismic and hydrologic inverse problems. We found that the seismic inverse problem is better suited to the quantum annealer than the hydrologic inverse problem. This is because the ratio between a fast P-wave velocity and a slow P-wave velocity is small compared to the ratio between a high hydraulic conductivity and a low hydraulic conductivity. This ratio ultimately influences the variability of the coefficients in the Hamiltonian used to program the quantum annealer, with a large ratio resulting in higher variability. High variability in the Hamiltonian coefficients leads to poor performance because the small coefficients effectively get lost in the noise.

In addition to identifying a subsurface problem that is well-suited to the quantum annealer, we also developed methods that enable the quantum annealer to solve inverse problems with a realistic set of observations. This is a significant step forward because previous work, which focused on the hydrologic inverse problem, was limited to an unrealistic set of observations. In particular, it required that the hydraulic head be observed at every point on the computational grid. This was consistent with early methods that were used in computational hydrology—called direct inverse methods (Yeh,

By transitioning from hydrology to seismology and from a direct inverse method to an iterative inverse method, we have taken two significant steps toward enabling the use of quantum annealing for practical applications in subsurface geoscience. One significant hurdle remains, and that is increasing the resolution of the subsurface image that the quantum annealer can handle. This would be aided by adding additional qubits to the quantum annealer and increasing the qubits' connectivity, both anticipated in D-Wave's next quantum annealer.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: Greer (

DO'M conceived of the presented idea, with contributions from SG. SG implemented the method and performed the computations with support from DO'M. All authors discussed the results and contributed to the final manuscript. All authors contributed to the article and approved the submitted version.

SG acknowledges support from the United States Department of Energy through the Computational Science Graduate Fellowship (DOE CSGF) under grant number DE-SC0019323. DO'M acknowledges support from the National Nuclear Security Administration's Advanced Simulation and Computing program.

We would like to acknowledge John Golden for useful discussions.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.