Early steps toward practical subsurface computations with quantum computing

1. Introduction

Computation has played a critical role in subsurface geoscience for decades. It has been used to simulate ocean circulation (Pinardi et al., 1997), flow in subsurface fractures (Kosakowski and Berkowitz, 1999), mantle flow in stunning detail (Stadler et al., 2010), and has been instrumental to the modern developments in seismic imaging (Bednar, 2005), among many other examples. Computers are also widely used to estimate subsurface properties that are difficult to observe directly using inverse analysis (Khan et al., 2000; Lu and Robinson, 2006). Massive performance improvements have buoyed the widespread use of computers.

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, 2018) and machine learning (Gentine et al., 2018), either together or separately. Approaches such as these have the potential to improve performance significantly, but only by a constant factor. That is, they might provide a speed-up of a factor of, e.g., 10 or 100, but the speed-up does not depend on the problem size. Quantum computing, on the other hand, opens the doors to fundamentally different algorithms that give larger and larger speed-ups as the problem becomes larger. This improved scaling behavior is crucial in subsurface geoscience. Solving equations, such as the groundwater flow equation or the seismic wave equation, might require a highly-refined mesh to resolve heterogeneity. This is because heterogeneities can be very small.

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., 2019; Zhong et al., 2020; Morvan et al., 2023). While it is known that some problems related to cryptography and quantum chemistry are well-suited to quantum computers, the picture is much less clear for computational science broadly. There are efficient algorithms that could theoretically be used to solve large systems of equations (Harrow et al., 2009; Subaşı et al., 2019), but these efficiencies may be undone by the implementation details needed to use these algorithms for a particular application (Aaronson, 2015). Current work in the geosciences utilizes a more empirical approach—trying different problems and observing the performance (O'Malley, 2018; Sarkar and Levin, 2018; Greer and O'Malley, 2020; Dukalski, 2021; Henderson et al., 2021; Souza et al., 2022; Dukalski et al., 2023). Generally, the performance on current quantum computers lags behind the performance of classical computers using the best algorithms. Nonetheless, it is a critical first step to establish whether or not a quantum computer can solve a problem before trying to establish a performance advantage.

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., 2022). Past work in hydrology required the use of an unrealistic set of observations (O'Malley, 2018), whereas the approach used here can handle arbitrary, realistic sets of observations. While we study this approach in the context of hydrologic and seismic inverse analysis, it can be applied to other subsurface applications where the goal of the inverse analysis is to segment the subsurface into two separate facies. Traditionally, a seismic inverse problem of this nature could be solved using imaging methods such as reverse time migration (Baysal et al., 1983) or full waveform inversion (Virieux and Operto, 2009). Hydrologic inverse problems are usually solved using variants of the geostatistical approach such as the principal component geostatistical approach (Kitanidis and Lee, 2014) and sometimes more modern techniques leveraging machine learning are used (Kadeethum et al., 2021; Wu et al., 2023).

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., 2012). Some quantum algorithms show an exponential speed-up (e.g., Harrow et al., 2009; Wiebe et al., 2012) which could be transformational in this context. A quadratic speedup [often built upon Grover's algorithm (Grover, 1996)] is more common and could still be impactful for subsurface inverse analysis, which is often computationally expensive. These algorithms could open doors to solve subsurface geoscience problems with unprecedented resolution and accuracy. After decades of research, the improvements that can be made in classical computational methods are becoming marginal (Shalf, 2020). Novel computational architectures, of which quantum computing is arguably the most promising, remain largely unexplored and have tremendous potential. Now is the time to do this exploration.

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.

2. Methods

2.1. Quantum annealing

Quantum annealing is a heuristic optimization algorithm, similar to simulated annealing (Kirkpatrick et al., 1983), that seeks to find optimal solutions faster than classical methods by exploiting quantum fluctuations (Kadowaki and Nishimori, 1998). There are theoretical guarantees of convergence to the optimal state under certain conditions (Morita and Nishimori, 2008). In practice, these assumptions are generally violated. For example, with the D-Wave quantum annealers (Johnson et al., 2011) that we use here, the anneal process is often performed quickly, whereas the theoretical guarantees generally require the anneal to be performed slowly. In the language of quantum annealing, this means that the adiabaticity is violated, since the fast annealing process means the system will often leave the ground state.

The input to a D-Wave quantum annealer is a vector h = (h_i) and a matrix J = (J_ij). The matrix, J, is sparse with a sparsity pattern defined by the connectivity graph associated with the annealer's qubits. The size of the matrix, N, is determined by the size of the problem. For existing D-Wave quantum annealers, this is based on a so-called Chimera graph (Boothby et al., 2016).

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

$\begin{array}{ll}g(s)={\displaystyle \sum _{i=1}^N}{h}_i{s}_i+{\displaystyle \sum _{i,j=1}^N}{J}_{ij}{s}_i{s}_j,& (1)\end{array}$

where each spin, s_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 J_ij. That is, some of the J_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 h_i and J_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 s that make f(s) small. This distribution can often be well-approximated by a Boltzmann distribution where Equation 1 defines the energy subject to a temperature or “effective temperature” that is a characteristic of the hardware. Equation 1 can be formulated as a quadratic unconstrained binary optimization (QUBO) problem,

$\begin{array}{ll}g(q)={\displaystyle \sum _{i=1}^N}{a}_i{q}_i+{\displaystyle \sum _{i,j=1}^N}{b}_{ij}{q}_i{q}_j& (2)\end{array}$

where each bit, q_i, is either 0 or 1, and is related to the formulation in Equation 1 via s_i = 2q_i − 1. This change, from inverting values of s_i ∈ {−1, 1} as in Equation 1 to q_i ∈ {0, 1} as in Equation 2, is what differentiates the Ising model from the QUBO model. The values of h_i and J_ij can also be transformed into the values of a_i and b_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.

2.2. Inverse approach

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

$\begin{array}{ll}F(c)={\displaystyle \sum _{i=1}^N}|{\left[U(c)\right]}_i-{û}_i{|}^2,& (3)\end{array}$

where U is a non-linear forward modeling operator such that U(c) = u is a solution to the relevant governing equation (wave equation for the seismic problem, groundwater flow equation for the hydrology problem), c is the subsurface model, $\widehat{u}$ is the measurements (wavefield in the case of the seismic problem and hydraulic head in the hydrology problem), and i represents the index of the N different observations. Note that each component of c can take only two values, either c_low or c_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, c⁽⁰⁾, and iteratively produce c^(k+1) from c^(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

$\begin{array}{ll}F(q)={\displaystyle \sum _{n=1}^N}|{\left[\tilde{U}(q)\right]}_n-{\widehat{u}}_n{|}^2& (4)\end{array}$

where $\tilde{U}(q)$ is a linear approximation to U at the current best estimate c^(k). This estimate is given by

$\begin{array}{ll}\tilde{U}{(q)}^{(k)}=U(c^{(k)})+{\displaystyle \sum _{m=1}^M}{q}_m{B}_m& (5)\end{array}$

where q is a binary vector that indicates how the model should be updated, and M is the size of the model, which is the number of parameters that we estimate. The vector B_m is defined to be

$\begin{array}{ll}{B}_m=\{\begin{array}{ll}U(c^{(k)}+[c_{low}-c_{high}]e_i)-U(c^{(k)})\hfill & \text{if }{c}_i^{(k)}=c_{high}\hfill \\ U(c^{(k)}+[c_{high}-c_{low}]e_i)-U(c^{(k)})\hfill & \text{if }{c}_i^{(k)}=c_{low}\hfill \end{array}& (6)\end{array}$

where e_i is the i^th standard basis vector of size M consisting of all zeros except for a 1 in the i^th component. Essentially, B_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 M forward model runs. When updating the model, if q_i = 0 then $c_i^{(k+1)}$ takes the same value as $c_i^{(k)}$. On the other hand, if q_i = 1 then $c_i^{(k+1)}$ takes the same opposite value as $c_i^{(k)}$. That is, if q_i = 1 then

$\begin{array}{ll}{c}_i^{(k+1)}=\{\begin{array}{ll}{c}_{high}\hfill & \text{if }{c}_i^{(k)}=c_{high}\text{ and }{q}_i=0\hfill \\ c_{low}\hfill & \text{if }{c}_i^{(k)}=c_{low}\text{ and }{q}_i=0\hfill \\ c_{low}\hfill & \text{if }{c}_i^{(k)}=c_{high}\text{ and }{q}_i=1\hfill \\ c_{high}\hfill & \text{if }{c}_i^{(k)}=c_{low}\text{ and }{q}_i=1\hfill \end{array}& (7)\end{array}$

The least-squares objective function for iteration k is then equivalent to

$\begin{array}{ll}F(q)=\left|\right|U(c^{(k)})-\widehat{u}|{|}_2^2+2{(U(c^{(k)})-\widehat{u})}^T{B}^{T}q+q^{T}B{B}^{T}q.& (8)\end{array}$

Note that this equation has a constant term, a term that is linear in q, and a term that is quadratic in q. The constant term is irrelevant to the optimization process, and neglecting it results in a function of the form in Equation (2). After creating the QUBO for a given iteration (Equation 8), Los Alamos National Laboratory's D-Wave 2000Q quantum annealer does the (forward) annealing and returns 1,000 possible solutions. We analyze the first several solutions that minimize the local objective function (Equation 8), select the update among those that minimizes the global objective function (3), and update the model accordingly. The selected update, q ∈ {0, 1}^M, is used to update the model, c^(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.

2.2.1. Hydrologic inverse analysis

We study the steady-state groundwater flow equation,

$\begin{array}{ll}\nabla ·\left(k(x)\nabla h\right)=f& (9)\end{array}$

where h represents the hydraulic head, k(x) denotes the heterogeneous hydraulic conductivity, and f represents fluid sources or sinks. Note that throughout, we will assume that f = 0 (i.e., there are no fluid sources/sinks) and that k(x) is either k_h or k_l where k_h is a high conductivity value and k_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, k(x). That is, to determine at each location, x, whether k(x) = k_h or k(x) = k_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 k_h or k_l. We first look at the example from Lu and Robinson (2006), which includes two low-permeability zoned embedded in a high permeability background medium. Our second hydrology example includes a slurry wall in our domain of interest. The hydraulic conductivities are distributed on a grid that is coarse compared to the finite volume grid on which the Equation (9) is solved. This use of these two grids reduces the number of hydraulic conductivity variables so that they can all be fit on the quantum annealer and the finite volume grid remains sufficiently resolved to produce a physically accurate simulation. Given these hydraulic conductivities, Equation (9) is solved to obtain a set of hydraulic head observations. These observations and the current estimate of the hydraulic conductivities are then used to formulate a QUBO using the previously discussed approach. Effectively, the PDE solver used to solve Equation (9) provides the function F(c) in Equation (3). After approximations, this results in the QUBO given in Equation (8). The quantum annealer is then used to optimize the QUBO, providing an updated estimate of the hydraulic conductivities. This iterative process continues until the convergence criteria are satisfied.

We have a 700 × 700 meter domain in this application, with seven discrete permeability blocks in the x-direction and seven in the y-direction. We place 24 receivers across the surface of the domain in a checkerboard pattern. Our computational mesh grid is dx = dy = 10 meters. Our goal is to invert for the locations of the two different facies, which each have permeabilities of k_l or k_h. Given this problem geometry, there are 2⁴⁹ ≈ 5.6 × 10¹⁴ possible solutions. The exact model we use is shown in Figure 1.

FIGURE 1

Figure 1. 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.

2.2.2. Seismic inverse analysis

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,

$\begin{array}{ll}\left(\frac{1}{c{(x)}^2}\frac{{\partial }^2}{\partial t^2}-\nabla \right)u=f,& (10)\end{array}$

where c(x) is the P-wave velocity, u is the measured wavefield, and f is the forcing, or source term. Similar to the hydrology example, we assume that c(x) is either c_high or c_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 f to be a 25 Hz Ricker wavelet at the top center of the domain. We also spread seven receivers evenly across the surface of the domain to record the wavefield measurements. We use a 1 millisecond sampling interval and record for 0.4 s at all sensor locations for the wavefield measurements. We use the observations recorded at the receiver locations to formulate a QUBO using the method discussed in this paper. As in the hydrology problem, the PDE solver used to solve Equation (10) provides the function F(c) in Equation (3). After approximations, this results in the QUBO given in Equation (8). The exact velocity model we use is in Figure 2. Given this problem geometry, there are 2⁵⁰ ≈ 1.1 × 10¹⁵ possible solutions.

FIGURE 2

Figure 2. 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 x-direction and five in the z-direction. The domain is 1, 000 m wide and 500 m deep, so each discrete velocity block is 100 × 100 m. The computational mesh grid is dx = dz = 10 meters. We have seven receivers, denoted by red triangles, which are spread evenly across the top of the domain. The source is located at the top center of the domain.

3. Results

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 c_high/c_low is large and another where c_high/c_low is relatively small. For the seismic problem, the case where c_high/c_low is relatively small is realistic. On the other hand, the case where c_high/c_low is large is realistic for the hydrologic problem.

Figure 3 shows the convergence behavior of the inverse analyses for the hydrologic problem. The left panel shows the convergence pattern when there is low contrast between the high hydraulic conductivity and low conductivity, while the right panel shows the convergence pattern when there is a large contrast. For the low contrast case, we use $k_{low}=1\times 10^{-3}m/s$ and $k_{high}=2\times 10^{-3}m/s$, and for the high contrast case, we use $k_{low}=5\times 10^{-8}m/s$ and $k_{high}=5\times 10^{-3}m/s$. In both cases, the initial model c₀ = k_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.

FIGURE 3

Figure 3. The convergence of our quantum annealing inverse approach for 10 low (left) and high (right) contrast hydrologic inverse problems is shown using the same initial value of c₀ = k_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.

Figure 4 shows the convergence behavior of the inverse analyses for the seismic problem. The left panel shows the convergence pattern when there is a low contrast between the high velocity and low velocity, while the right panel shows the convergence pattern when there is a large contrast.

FIGURE 4

Figure 4. The convergence of our quantum annealing inverse approach for 10 low (left) and high (right) contrast seismic inverse problems is shown using the same initial value of c₀ = c_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: c_low = 4, 250m/s and c_high = 4, 750m/s. In the high contrast case, we use c_low = 2, 000m/s and c_high = 5, 000m/s. In both cases, the initial model c₀ = c_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.

4. Discussion

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. Table 1 shows hydraulic conductivity ranges for a set of unconsolidated sediments (Fetter, 2018). Note that even within one class of materials, the hydraulic conductivity can vary by several orders of magnitude. Since different materials often coexist in the same region of the subsurface, the variability can be even larger than this. For example, the widely used SPE 10 model exhibits variation in the permeability over 8–12 orders of magnitude from 10⁻⁷ milliDarcy to 10⁴ milliDarcy (Lie, 2019).

TABLE 1

Table 1. Variability in permeability in unconsolidated sediments (Fetter, 2018).

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. Figure 5 shows the distribution of the coefficients in one iteration of the seismic and hydrologic inverse analysis for both the low and high contrast cases. The coefficients that are small in magnitude tend to be lost in the noise associated with the quantum annealer (Golden and O'Malley, 2021). This is caused by the rescaling that is necessary to fit the coefficients within the dynamic range allowed by the quantum annealer. The quantum annealer has little or no information about qubits whose linear and quadratic coefficients are below the hardware's analog noise level. Since the iterative model updates are determined by solving this QUBO problem, which may be inaccurate in high contrast cases due to its coefficients being below the threshold of hardware noise, the model updates selected from this method may diverge from minimizing the original problem's objective function, as seen in Figure 3.

FIGURE 5

Figure 5. The cumulative distribution function for an Ising problem in (A) the seismic inverse analysis and (B) the hydrologic inverse analysis are shown. Note that the high contrast problem is more likely to have small coefficients that get lost in the noise of the quantum annealer.

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 Table 2. In general, the P-wave velocity of a given unit will increase with depth since the material becomes more compact. Because of this, units close to each other in depth are likely to have a more similar velocity than units far from each other in depth. This allows for a lower impedance contrast under the assumption of a constant-density acoustic model. The high-contrast case, where the performance is not as good, would be uncommon in real examples of subsurface velocity models with the notable exception in areas with salt bodies, such as in the Gulf of Mexico. Because of the difficulties with QUBO coefficients being lost in the noise, quantum annealing appears more suitable for seismic inverse analysis than hydrologic inverse analysis using the method we propose.

TABLE 2

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

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., 2019). This hardware significantly increases the number of qubits on the quantum annealer chip and the number of connections per qubit. The effect of this will be to approximately triple the number of parameters that can be calibrated. Approximately tripling the number of parameters will significantly increase the resolution of the inverse model that the quantum annealer can handle. The Advantage hardware also reduces the noise on the system, which could improve the performance of high contrast problems.

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.

5. Conclusion

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, 1986). The transition to an iterative approach for inverse analysis with quantum annealing brings it in line with modern methods for inverse analysis that also use iterative methods and can handle realistic observation sets.

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.

Data availability statement

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 (2020).

Author contributions

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.

Funding

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.

Acknowledgments

We would like to acknowledge John Golden for useful discussions.

Conflict of interest

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.

Publisher's note

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.

References

Aaronson, S. (2015). Read the fine print. Nat. Phys. 11, 291–293. doi: 10.1038/nphys3272

CrossRef Full Text | Google Scholar

Arute, F., Arya, K., Babbush, R., Bacon, D., Bardin, J. C., Barends, R., et al. (2019). Quantum supremacy using a programmable superconducting processor. Nature 574, 505–510. doi: 10.1038/s41586-019-1666-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Baysal, E., Kosloff, D. D., and Sherwood, J. W. (1983). Reverse time migration. Geophysics 48, 1514–1524. doi: 10.1190/1.1441434

CrossRef Full Text | Google Scholar

Bednar, J. B. (2005). A brief history of seismic migration. Geophysics 70, 3MJ–20MJ. doi: 10.1190/1.1926579

CrossRef Full Text | Google Scholar

Boothby, T., King, A. D., and Roy, A. (2016). Fast clique minor generation in chimera qubit connectivity graphs. Quant. Inf. Process. 15, 495–508. doi: 10.1007/s11128-015-1150-6

CrossRef Full Text | Google Scholar

Bourbié, T., Coussy, O., and Zinszner, B. (1987). Acoustics of Porous Media. Editions TECHNIP.

Google Scholar

Dattani, N.S., Szalay, S., and Chancellor, N. (2019). Pegasus: The second connectivity graph for large-scale quantum annealing hardware. arXiv [Preprint]. arXiv: 1901.07636.

Google Scholar

Dukalski, M. (2021). “Toward an application of quantum computing in geophysics,” in Fifth EAGE Workshop on High Performance Computing for Upstream, Vol. 2021 (European Association of Geoscientists & Engineers), 1–5.

Google Scholar

Dukalski, M., Rovetta, D., van der Linde, S., Möller, M., Neumann, N., and Phillipson, F. (2023). Quantum computer-assisted global optimization in geophysics illustrated with stack-power maximization for refraction residual statics estimation. Geophysics 88, V75–V91. doi: 10.1190/geo2022-0253.1

CrossRef Full Text | Google Scholar

Fatemi, S., and Poppe, A. R. (2018). Solar wind plasma interaction with asteroid 16 psyche: implication for formation theories. Geophys. Res. Lett. 45, 39–48. doi: 10.1002/2017GL073980

CrossRef Full Text | Google Scholar

Fetter, C. W. (2018). Applied Hydrogeology. Long Grove, IL: Waveland Press.

Google Scholar

Gentine, P., Pritchard, M., Rasp, S., Reinaudi, G., and Yacalis, G. (2018). Could machine learning break the convection parameterization deadlock? Geophys. Res. Lett. 45, 5742–5751. doi: 10.1029/2018GL078202

CrossRef Full Text | Google Scholar

Golden, J. K., and O'Malley, D. (2021). Pre- and post-processing in quantum-computational hydrologic inverse analysis. Quantum Inf. Process 20, 176. doi: 10.1007/s11128-021-03115-y

CrossRef Full Text | Google Scholar

Greer, S. (2020). sygreer/QuantumAnnealingInversion.jl: First release of QuantumAnnealingInversion.jl (v1.0.0). Zenodo. Available online at: https://zenodo.org/record/4313834

Google Scholar

Greer, S., and O'Malley, D. (2020). “An approach to seismic inversion with quantum annealing,” in SEG Technical Program Expanded Abstracts (SEG), 2845–2849. doi: 10.1190/segam2020-3424413.1

CrossRef Full Text | Google Scholar

Grover, L. K. (1996). “A fast quantum mechanical algorithm for database search,” in Proceedings of the Twenty-Eighth Annual ACM Symposium on Theory of Computing (Association for Computing Machinery), 212–219.

Google Scholar

Harrow, A. W., Hassidim, A., and Lloyd, S. (2009). Quantum algorithm for linear systems of equations. Phys. Rev. Lett. 103, 150502. doi: 10.1103/PhysRevLett.103.150502

PubMed Abstract | CrossRef Full Text | Google Scholar

Henderson, J. M., O'Malley, D., and Viswanathan, H. S. (2021). “Interrogating the performance of quantum annealing for the solution of steady-state subsurface flow,” in 2021 IEEE High Performance Extreme Computing Conference (HPEC) (IEEE), 1–6.

Google Scholar

Johnson, M. W., Amin, M. H., Gildert, S., Lanting, T., Hamze, F., Dickson, N., et al. (2011). Quantum annealing with manufactured spins. Nature 473, 194–198. doi: 10.1038/nature10012

PubMed Abstract | CrossRef Full Text | Google Scholar

Kadeethum, T., O'Malley, D., Fuhg, J. N., Choi, Y., Lee, J., Viswanathan, H. S., et al. (2021). A framework for data-driven solution and parameter estimation of pdes using conditional generative adversarial networks. Nat. Comp. Sci. 1, 819–829. doi: 10.1038/s43588-021-00171-3

CrossRef Full Text | Google Scholar

Kadowaki, T., and Nishimori, H. (1998). Quantum annealing in the transverse ising model. Phys. Rev. E 58, 5355. doi: 10.1103/PhysRevE.58.5355

PubMed Abstract | CrossRef Full Text | Google Scholar

Khan, A., Mosegaard, K., and Rasmussen, K. L. (2000). A new seismic velocity model for the moon from a monte carlo inversion of the apollo lunar seismic data. Geophys. Res. Lett. 27, 1591–1594. doi: 10.1029/1999GL008452

CrossRef Full Text | Google Scholar

Kirkpatrick, S., Gelatt, C. D., and Vecchi, M. P. (1983). Optimization by simulated annealing. Science 220, 671–680. doi: 10.1126/science.220.4598.671

PubMed Abstract | CrossRef Full Text | Google Scholar

Kitanidis, P. K., and Lee, J. (2014). Principal component geostatistical approach for large-dimensional inverse problems. Water Resour. Res. 50, 5428–5443. doi: 10.1002/2013WR014630

PubMed Abstract | CrossRef Full Text | Google Scholar

Kosakowski, G., and Berkowitz, B. (1999). Flow pattern variability in natural fracture intersections. Geophys. Res. Lett. 26, 1765–1768. doi: 10.1029/1999GL900344

PubMed Abstract | CrossRef Full Text | Google Scholar

Lie, K. A. (2019). An Introduction to Reservoir Simulation Using MATLAB/GNU Octave: User Guide for the MATLAB Reservoir Simulation Toolbox (MRST). Cambridge: Cambridge University Press.

Google Scholar

Lu, Z., and Robinson, B. A. (2006). Parameter identification using the level set method. Geophys. Res. Lett. 33, L06404. doi: 10.1029/2005GL025541.1

CrossRef Full Text | Google Scholar

Morita, S., and Nishimori, H. (2008). Mathematical foundation of quantum annealing. J. Math. Phys. 49, 125210. doi: 10.1063/1.2995837

CrossRef Full Text | Google Scholar

Morvan, A., Villalonga, B., Mi, X., Mandrá, S., Bengtsson, A., Klimov, P. V., et al. (2023). Phase transition in random circuit sampling. arXiv [Preprint]. arXiv: 2304.11119. doi: 10.48550/arXiv.2304.11119

CrossRef Full Text | Google Scholar

O'Malley, D. (2018). An approach to quantum-computational hydrologic inverse analysis. Sci. Rep. 8, 1–9. doi: 10.1038/s41598-018-25206-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Pinardi, N., Korres, G., Lascaratos, A., Roussenov, V., and Stanev, E. (1997). Numerical simulation of the interannual variability of the mediterranean sea upper ocean circulation. Geophys. Res. Lett. 24, 425–428. doi: 10.1029/96GL03952

CrossRef Full Text | Google Scholar

Sarkar, R., and Levin, S. A. (2018). “Snell tomography for net-to-gross estimation using quantum annealing,” in SEG International Exposition and Annual Meeting (SEG).

Google Scholar

Shalf, J. (2020). The future of computing beyond moore's law. Philos. Transact. R. Soc. 378, 20190061. doi: 10.1098/rsta.2019.0061

PubMed Abstract | CrossRef Full Text | Google Scholar

Souza, A. M., Martins, E. O., Roditi, I., Sá, N., Sarthour, R. S., and Oliveira, I. S. (2022). An application of quantum annealing computing to seismic inversion. Front. Phys. 9, 748285. doi: 10.3389/fphy.2021.748285

CrossRef Full Text | Google Scholar

Stadler, G., Gurnis, M., Burstedde, C., Wilcox, L. C., Alisic, L., and Ghattas, O. (2010). The dynamics of plate tectonics and mantle flow: from local to global scales. Science 329, 1033–1038. doi: 10.1126/science.1191223

PubMed Abstract | CrossRef Full Text | Google Scholar

Subaçı, Y., Somma, R. D., and Orsucci, D. (2019). Quantum algorithms for systems of linear equations inspired by adiabatic quantum computing. Phys. Rev. Lett. 122, 060504. doi: 10.1103/PhysRevLett.122.060504

PubMed Abstract | CrossRef Full Text | Google Scholar

Virieux, J., and Operto, S. (2009). An overview of full-waveform inversion in exploration geophysics. Geophysics 74, WCC1–WCC26. doi: 10.1190/1.3238367

CrossRef Full Text | Google Scholar

Wiebe, N., Braun, D., and Lloyd, S. (2012). Quantum algorithm for data fitting. Phys. Rev. Lett. 109, 050505. doi: 10.1103/PhysRevLett.109.050505

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, H., Greer, S. Y., and O'Malley, D. (2023). Physics-embedded inverse analysis with algorithmic differentiation for the earth?s subsurface. Sci. Rep. 13, 718. doi: 10.1038/s41598-022-26898-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Yeh, W. W.-G. (1986). Review of parameter identification procedures in groundwater hydrology: the inverse problem. Water Resour. Res. 22, 95–108. doi: 10.1029/WR022i002p00095

CrossRef Full Text | Google Scholar

Zhong, H.-S., Wang, H., Deng, Y.-H., Chen, M.-C., Peng, L.-C., Luo, Y.-H., et al. (2020). Quantum computational advantage using photons. Science 370, 1460–1463. doi: 10.1126/science.abe8770

PubMed Abstract | CrossRef Full Text | Google Scholar