Application of Riemannian Seismic Ray Path Tracing in Salt Dome Prospecting

Abstract

Within the geophysical exploration utilising seismic methods, it is well known that if the explored distances are much greater than the wavelength of the seismic waves with which the exploration is carried out, the ray approach of the wave theory can be used. In this way, when the rays travel through an inhomogeneous medium, they follow curved trajectories, which is imperative to determine the geological features that produce reflection and refraction phenomena. In this paper, a simple algorithm for the calculation of the trajectory of a seismic beam through an inhomogeneous stratum is presented. For this, the construction of a pseudo-Riemannian metric is required from the function of P-wave velocities of the geological stratum. Thus, the problem is inverted because instead of finding the curved trajectory of the seismic beam in a background with a Euclidean metric, it is proposed that the beam follows a geodesic of a curved space-time specific to each stratum, becoming a simple and automatic process using the differential geometry apparatus. For the reader to gain insight into this tool, different geological setups from idealised ones up to a salt dome are presented.

1. Introduction

Ray theory is one of the most widely used methods in seismology and seismic exploration for the forward and inverse modelling of high-frequency seismic body waves. It provides a physical view of the wave propagation process by separating the wave field into individual elementary waves. This allows their identification [1] and path tracing of energy propagation through a specific medium [2]. However, it has some limitations, since it is approximate.

To elaborate further on the limitations of ray theory, it is crucial to consider its foundational assumptions and their implications. Firstly, ray theory assumes that seismic waves propagate as narrow beams along distinct paths, neglecting wave interactions and diffraction effects that occur in real-world geological environments. This simplification can lead to inaccuracies in predicting wave behaviour, especially in heterogeneous subsurface conditions where wave paths can be complex and irregular [3,4].

Secondly, ray theory often struggles with accurately modelling wave phenomena at lower frequencies or in regions with strong velocity gradients. These situations challenge the assumption of straight-line ray paths and can result in significant errors in seismic imaging and interpretation [3,4].

Moreover, while ray theory provides valuable insights into wave propagation dynamics, it does not account for wave phenomena such as scattering and attenuation, which are crucial for a more detailed understanding of seismic data. Also, it does not examine the amount of scattered energy from ray propagation throughout surfaces [5]. These omissions can limit the theory’s applicability in practical seismic exploration scenarios where comprehensive and precise modelling is essential [3,4].

While ray theory remains a cornerstone in seismic research and exploration, its approximations and inherent limitations necessitate careful consideration and often, supplementary methods to achieve accurate seismic imaging and interpretation in diverse geological settings.

Starting from an intuitive description of seismic propagation along special trajectories (rays), based on the laws of geometric optics [6], it has become a highly sophisticated method. The first seismological applications date back to the end of the 19th century, through the application of kinematic laws and straight rays to estimate travel times [3]. Jeffreys [7] was one of the pioneers in developing dynamical ray theory, also focusing on amplitudes and waveforms. The “modern” series theory (of rays) is obtained by calculating solutions to the equations of elastodynamics [8], and were obtained for the first time by Babich and Alekseyev [9] and afterwards by Karal Jr and Keller [10] for non-homogeneous isotropic media, and later by Babich [11] for anisotropic media.

Babich [11] proposes that high-frequency waves propagate with an oscillatory time dependence, which induces a solution of the wave equation by separation of variables. For the temporal part, a function for the travel time is made manifest, and the spatial part is expanded into a series of inverse powers of the angular frequency, whose coefficients are vector amplitude coefficients and are usually complex. For perfectly elastic media, travel times and vector amplitude coefficients are usually independent of frequency [3]. By introducing series expansion to the equations of elastodynamics, approximations of different orders to the ray method are obtained. However, the complexity of the equations increases considerably with increasing order of approximation, since higher-order amplitude coefficients are inaccurate and unstable, due to their high sensitivity to the fine details of the medium [3]. In this context, sensitivity refers to the extent to which higher-order amplitude coefficients depend on small variations in the properties of the medium. This means that even very small changes or imperfections in the medium can cause large variations in the calculated amplitude coefficients. As a result, higher-order approximations become less accurate and more prone to errors, making them unstable and difficult to use reliably in practical applications. In this way, only the zero-order coefficient, and at most the first-order coefficient, have been used in seismological applications. This fact has limited the use of the ray approximation in its conventional form. There are other seismic ray series such as Chapman [12], which expand in terms of particle velocity and traction, but they have not been popularised for ray tracing.

The zero-order ray method allows us to arrive at the well-known Eikonal equation, a nonlinear first-order partial differential equation for the travel time, the polarisation vector (which indicates the direction of the displacement vector) [3], and the transport equation that allows calculating amplitudes [13,14]. The Eikonal equation can be expressed by the Hamiltonian formalism [3,15], allowing it to be decomposed into a series of first-order nonlinear ordinary differential equations, often simplifying its (numerical) solution with appropriate initial conditions, obtaining both the trajectories and the travel times [16]. It is clear that by applying the corresponding Legendre transforms, one can arrive at a Lagrangian formulation of the seismic ray theory [1]. Furthermore, it is possible to arrive at these formulations through variational principles, such as those of Lagrange and Hamilton [17].

From this formalism, there are very different approaches such as paraxial ray theory [15,18], the study of the effects of structural interfaces [16,19,20], the study of the Green’s function in inhomogeneous and weakly anisotropic elastodynamics [21], perturbation methods in ray theory [4,12,22,23,24,25,26,27], handling around certain singularities [28,29,30,31,32], and even the study of “chaotic” rays using Lyapunov exponents [4,33]. It should be noted that ray theory is usually not applicable to S-type shear waves, which require separate treatments [26,34,35].

Currently, the ray method can be used only for high-frequency signals and in models where heterogeneity does not exceed a certain degree. Thus, it is desirable to have available a more general extension of ray methods, both to simplify the theory and to lead to more precise results. In this sense, there is Maslov’s asymptotic ray theory introduced to seismology by Chapman [12], Chapman and Drummond [36], Thomson and Chapman [37], the theory based on the sum of Gaussian rays [38,39], Maslov’s theory and Gaussian rays [40], the Gaussian packet sum method [4], the Kirchhoff surface integral method [12], and the one-way wave equation approach [41]. However, these extensions of ray theory only circumvent some of its original problems, hence, a new approach is required.

Considering that conventional ray theory can be written using variational principles, one may wonder whether Fermat’s principle in seismic can be completely rewritten into mathematical language without any degree of approximation. In this sense, there is a transparent and physically direct formulation of geometrical optics that allows the optical Fermat principle to be translated into the formalism of differential geometry, an area of mathematics that provides intrinsically curved spaces (manifolds) with the power of differential and integral calculus [42].

Carretero López et al. [43] show that the trajectories arising from Fermat’s principle behind geometrical optics are, mathematically speaking, equivalent to the geometrical problem of determining geodesic equations in an “optical space” that is not Euclidean (curved). In other words, the physical problem of determining the trajectories of a ray of light in the three-dimensional Euclidean space (${\mathbb{R}}^3$) is changed by the paradigm of finding the geodesics (minimum length trajectories, generalisation of the concept of line in ${\mathbb{R}}^3$) in an intrinsically curved mathematical space (differential manifold), which abstracts the essence of the material media in which these rays are propagated. In this way, all the tools of differential geometry are available to treat light rays, and therefore, this approach is constituted as a very robust mathematical approach to geometric optics.

Since electromagnetic and optical phenomena are essentially undulatory, the parallelism in the theory of rays is impeccable, finding that Fermat’s principle governs the phenomena of rays in both formalisms as they are identical in both cases. In this way, the demonstration by Carretero López et al. [43] that the trajectories arising from the Fermat’s principle are equivalent to the problem of determining the geodesic equations in a non-Euclidean (curved) “optical space” can be extended to seismics.

The purpose of this article is to extend the parallelism between geometrical optics and the differential geometry of Carretero López et al. [43] to the theory of seismic rays, where the main postulate will be that a seismic ray follows a geodesic curve (i.e., of minimum length) in an intrinsically curved “seismic space”. In equivalence to the optical refractive index, the curvature of such space will be determined by the specific inhomogeneities of the elastic propagation medium, for which purpose the inverse of the propagation speed of P waves in geological media is proposed as “refractive index of seismic rays”.

For this purpose, Section 2 presents the elements of differential geometry required to develop this work, while in Section 3.1 it is applied to construct the Riemannian seismic ray theory. In order to summarise, Section 3.2 presents the methodology to be followed to apply this theory, while Section 4 presents several relevant applications in geophysics, such as a stratum change at a certain depth, simplified modelling of the speed of P-wave in the mantle, as well as an example of the modelling of irregular interfaces between two strata of constant wave propagation velocity. These simple examples show that the theory presented here reproduces the expected results. It should be noted that the last example shows the application of the theory to a realistic model of the subsoil with speeds of great interest in geophysical exploration, which contains a salt dome. In this case, it can be clearly seen that the seismic rays are refracted in the salt dome due to their high contrast with respect to their surroundings, which is consistent with the shadow generated by these objects in exploration and geophysical simulation. Finally, Section 5 presents an extensive discussion of the results and some interesting concluding observations are provided in Section 6.

2. Materials and Methods (Theoretical Background)

Differential geometry is a profound and intricate area of mathematics that extends the concepts and tools of differential and integral calculus to spaces with intrinsic curvature. This field is essential for understanding and analysing the geometric properties of various mathematical and physical systems. It involves the study of curves, surfaces, and higher-dimensional manifolds, which are spaces that locally resemble Euclidean space but can have complex global structures. Below is a detailed exploration of key concepts in differential geometry. For a more rigorous and comprehensive presentation of these concepts, it is recommended to review the texts of [44,45,46,47,48,49,50,51,52,53,54,55]. For readers more interested in a geometrical approach, they can consult [44,46,49], while those seeking a more tensor-like approach can review [45,48,50,52]. Finally, for readers interested in a presentation more linked to physical phenomena, they can consult [47,51]. In addition [53] for a discussion more focused in mechanics and [54,55] for those more interested from an astrophysics/cosmology perspective.

2.1. Differential Manifold

A manifold is essentially a continuous space that locally resembles the Euclidean space. This space is a set that can be continuously parametrized, the number of independent parameters being the dimension of the manifold. Think of $\mathbb{M}$, a d-dimensional manifold, the d parameters that parametrize it, ${\left\{q^i\right\}}_{i=1}^d$, are global functions in $\mathbb{M}$, and are a coordinate system for said manifold. In the same way, a curve $\kappa \in \mathbb{M}$ is a smooth function of $\kappa :I\subset \mathbb{R}\to \mathbb{M}$. and given a curve $\kappa \in \mathbb{M}$, the vector tangent $\overrightarrow{U}$ to such curve at each point p of $\mathbb{M}$ can be defined; the set of these vectors forms the space tangent to $\mathbb{M}$ at the point p, $T_p\left(\mathbb{M}\right)$,which is a vector space. In the same way, a tensor field assigns to each point of the space a tensor M-contravariant and N-covariant, or type $\left(\begin{array}{c}M\\ N\end{array}\right)$.

A differential manifold is one that has a differentiable structure, i.e., a continuous and differentiable space. This implies that at each point of the manifold one can define a scalar field $\varphi$ that can be differentiated anywhere. In other words, in each curve $\kappa$ in $\mathbb{M}$ (parametrized by, say, $\lambda$) is defined $d\varphi /d\lambda ={\nabla }_{\overrightarrow{U}}\varphi$, where the last equality follows by the vector field tangent to the curve.

2.1.1. Metric or Riemannian Manifold

A metric or Riemannian manifold is one where a notion of distance between any two points of the manifold is also implemented. This can be done using the definition of a global function in $\mathbb{M}$ (tensor $\left(\begin{array}{c}0\\ 2\end{array}\right)$).

Definition 1.

A Riemannian metric tensor $\mathbf{g}:\mathbb{M}\times M\to \mathbb{R}$ is a type $\left(\begin{array}{c}0\\ 2\end{array}\right)$ tensor that satisfies the following properties [56]:

1.

 $\mathbf{g}(X,Y)\ge 0$

2.

 $\mathbf{g}(X,Y)=0\iff X=Y$

3.

 $\mathbf{g}(X,Y)=\mathbf{g}(Y,X)$

4.

 $\mathbf{g}(X,Z)\le \mathbf{g}(X,Y)+\mathbf{g}(Y,Z)$

$\forall X,Y,Z\in T_p\left(\mathbb{M}\right)$.

The components of the metric tensor $\mathbf{g}$ with respect to the coordinate system $\left\{q^i\right\}$ are

$$g_{ij}=\mathbf{g}\left({\mathsf{\partial }}_i,{\mathsf{\partial }}_j\right),$$

(1)

where ${\mathsf{\partial }}_i:=MM\mathsf{\partial }/\mathsf{\partial }{q}_i$. The metric is a symmetric tensor, i.e., $g_{ij}=g_{ji}$.

Definition 2.

A metric or Riemannian space $\left(\mathbb{M},\mathbf{g}\right)$ is a manifold equipped with a metric tensor $\mathbf{g}$.

Covariant Derivative

Now it is important to define the differentiation in a Riemannian manifold. Clearly, in order to calculate the change of a vector along a curve (different points of a manifold), it should be considered that at each point $p_j\in \mathbb{M}$, $T_{p_j}\left(\mathbb{M}\right)$ changes, and that change is reflected at the base of $T_{p_j}\left(\mathbb{M}\right)$ at each point. This change, due to the curvature, must be considered when calculating the derivative of a vector (or, in general, of any tensor) along an intrinsically curved manifold. The correction factor is always proportional to a set of coefficients, called Christoffel symbols, which in the coordinate system $\left\{q^i\right\}$ can be determined through the components of the metric tensor (1) as

$${\Gamma }_{mn}^l=\frac{1}{2}\sum _{k=1}^d{g}^{lk}\left({\mathsf{\partial }}_m{g}_{nk}+{\mathsf{\partial }}_n{g}_{mk}-{\mathsf{\partial }}_k{g}_{mn}\right),$$

(2)

where $l,m,n=1,\cdots ,d$, $g^{lk}$ will be the components of the inverse matrix of the metric coefficients (1). The Christoffel symbols are symmetrical in their two covariant indices, i.e., ${\Gamma }_{mn}^l={\Gamma }_{nm}^l$. In this way, the derivative of a vector will be given as follows.

Theorem 1.

Let $\kappa :I\subset \mathbb{R}\to \mathbb{M}$ be a curve in $\left(\mathbb{M},\mathbf{g}\right)$ (parametrized by λ), and let ${{\rm Y}}_{\kappa }\left(\mathbb{M}\right)$ be the space of all the smooth vector fields induced by κ. The covariant derivative along κ will be the function:

$${\displaystyle \frac{D}{D\lambda }}:{{\rm Y}}_{\kappa }\left(\mathbb{M}\right)\to {{\rm Y}}_{\kappa }\left(\mathbb{M}\right),$$

(3)

which satisfies the properties of linearity, additivity and the Leibniz rule as well.

This definition can be generalised trivially for the derivative of any tensor type $\left(\begin{array}{c}M\\ N\end{array}\right)$. In the coordinate system $\left\{q^i\right\}$, the derivative (3) of the components of the vector $\overrightarrow{V}$ vector will be given by

$$\frac{D{V}^i}{D\lambda }=\sum _{j=1}^d{U}^j{\mathsf{\partial }}_j{V}^i+\sum _{j,k=1}^d{\Gamma }_{kj}^i{V}^k{U}^j,$$

(4)

where $\overrightarrow{U}=\frac{d\overrightarrow{q}}{d\lambda }$ is the vector tangent to the curve $\kappa$.

2.1.2. Geodesic Equations

A geodesic is the generalisation of the straight line in a curved manifold, and is defined as:

Definition 3.

A curve $\kappa :I\subset \mathbb{R}\to \mathbb{M}$ is a geodesic curve in $\left(\mathbb{M}\right)$ if the covariant derivative of its parametric velocity is zero, i.e.,

$$\frac{D\dot{\kappa }}{D\lambda }=0.$$

(5)

But since $\dot{\kappa }:=\overrightarrow{U}=\frac{d\overrightarrow{q}}{d\lambda }$, in the coordinate system $\left\{q^i\right\}$, this condition is equivalent to

$$\frac{D{q}^i}{D\lambda }=\frac{d}{d\lambda }\left(\frac{d{q}^i}{d\lambda }\right)+\sum _{k=1}^d{\Gamma }_{mn}^k\frac{d{q}^m}{d\lambda }\frac{d{q}^n}{d\lambda }=0$$

(6)

and since the Christoffel symbols ${\Gamma }_{mn}^l$ are known functions of the $\left\{q^i\right\}$ coordinates, this is a nonlinear (quasi-linear) second-order differential equation for $q^i=q^i\left(\lambda \right)$ [57].

Although for this work it is sufficient and necessary to know the formalism of the geodesic equations, below are the tools of differential geometry that deal with the concepts of intrinsic curvature. These will allow us to gain better intuition of the effect of the curvature in the seismic rays.

2.2. Curvature and Curvature Tensor

To determine whether a space is intrinsically curved or not, it is necessary to evaluate the components of the Riemann curvature tensor, which characterises the curvature in a tensor way. Under a coordinate system $\left\{q^i\right\}$ in $\mathbb{M}$, the components of the Riemann $\left(\begin{array}{c}1\\ 3\end{array}\right)$ tensor take the form

$$R_{lmn}^k={\mathsf{\partial }}_m{\Gamma }_{ln}^k-{\mathsf{\partial }}_n{\Gamma }_{lm}^k+\sum _{i=1}^d{\Gamma }_{im}^k{\Gamma }_{ln}^i-\sum _{i=1}^d{\Gamma }_{in}^k{\Gamma }_{lm}^i.$$

(7)

In this way, a flat manifold is one in which $R_{lmn}^k=0\forall k,l,m,n=1,\cdots ,d$. The tensor $R_{lmn}^k$ has $\left(d^4-d^2\right)/12$ independent components, since it satisfies the properties of symmetry $R_{abcd}=R_{cdab}$ and anti-symmetry $R_{abcd}=-R_{bacd}=-R_{abdc}$ as well as the identities of Bianchi. Moreover, any other tensor that describes the intrinsic curvature of a manifold can be deduced from this tensor. Such is the case of the Ricci tensor,

$$R_{mn}=\sum _{k=1}^d{R}_{mkn}^k,$$

(8)

and of the Ricci scalar field

$$R=\sum _{m,n=1}^d{g}^{mn}{R}_{mn}.$$

(9)

The latter is the simplest curvature invariant in a Riemannian manifold. Being a scalar invariant, it assigns a real number to each point of $\mathbb{M}$ according to the geometry of that manifold in the vicinity of the point, whose value does not depend on the coordinate system that parametrizes $\mathbb{M}$.

3. Development

3.1. Riemannian Geophysical Modelling

For the geophysical case, the differential manifolds will be generally of dimension $d=3$ or $d=2$, i.e., ${\mathbb{M}}^3$ or ${\mathbb{M}}^2$, respectively. Let a stratum of propagation velocity $v=v\left(q^i\right)$.

Definition 4.

The seismic refractive index $n:{\mathbb{M}}^d\to \mathbb{R}$ is defined as a function of class $C_{{\mathbb{M}}^d}^{\infty }$ such that

$$n=n\left(q^i\right)=\frac{1}{v}$$

(10)

Lemma 1.

If $n=n\left(q^i\right):{\mathbb{M}}^d\to \mathbb{R}$ is a smooth function,

$$\mathbf{g}=n^2\delta ,$$

(11)

is the Riemannian metric induced by n in ${\mathbb{M}}^d$, where δ is the Euclidean metric (with components ${\delta }_{ij}$, the Kronecker delta) [43].

Corollary 1.

The Riemannian seismic metric (11) naturally parametrizes the differential manifold $\mathbb{M}$ in a Cartesian coordinate system, i.e., $\left\{q^i\right\}$ are Cartesian coordinates.

Definition 5.

A seismic space $\left(\mathbb{M},\mathbf{g};n\right)$ is defined as a Riemannian manifold $\left(\mathbb{M},\mathbf{g}\right)$, with $\mathbf{g}$ the metric induced by the refraction index n.

Theorem 2.

(Fermat’s Principle) A curve κ in a seismic space $\left(\mathbb{M},\mathbf{g};n\right)$ is a seismic ray if and only if κ is a geodesic in that space, i.e., if

$$\frac{D\dot{\kappa }}{D\lambda }=0.$$

(12)

The equivalence between Fermat’s Principle and the statement of Theorem 2 is demonstrated by Carretero López et al. [43]. For a seismic space $\left(\mathbb{M},\mathbf{g};n\right)$, the geodesic Equation (12) produce seismic rays which follow the laws of transmission/refraction.

Definition 6.

Let $\left(\mathbb{M},\mathbf{g};n\right)$ be a seismic space and let $R_{lmn}^k$ be its Riemann curvature tensor. The reverse seismic space, or seismic space of inverse curvature, is defined as the seismic space $\left({\mathbb{M}}^{\prime },{\mathbf{g}}^{\prime };n^{\prime }\right)$ whose curvature tensor $R^{\prime }$ is such that

$${R^{\prime }}_{lmn}^k=-{R^k}_{lmn}$$

(13)

$\forall k,l,m,n=1,\cdots ,d$.

Theorem 3.

Let $\left(\mathbb{M},\mathbf{g};n\right)$ be a seismic space of dimension $d=2$. The seismic space induced by the inverse of the refractive index, $n^{-1}$, i.e., $\left(\mathbb{M},{\mathbf{g}}^{-1};n^{-1}\right)$, is the seismic space of inverse curvature to $\left(\mathbb{M},\mathbf{g};n\right)$ of dimension $d=2$.

Proof. 

As $\left(\mathbb{M},\mathbf{g};n\right)$ is a seismic space, in a Cartesian coordinate system ${\left\{q^i\right\}}_{i=1}^2$, according to Equation (2), induces the following Christoffel symbols:

$${\Gamma }_{11}^1=\frac{{\mathsf{\partial }}_{1}n}{n}{\Gamma }_{12}^1={\Gamma }_{21}^1=\frac{{\mathsf{\partial }}_{2}n}{n}{\Gamma }_{22}^1=-\frac{{\mathsf{\partial }}_{1}n}{n}$$

$${\Gamma }_{11}^2=-\frac{{\mathsf{\partial }}_{2}n}{n}{\Gamma }_{12}^2={\Gamma }_{21}^2=\frac{{\mathsf{\partial }}_{1}n}{n}{\Gamma }_{22}^2=\frac{{\mathsf{\partial }}_{2}n}{n}$$

Naturally, because $d=2$, the curvature tensor (Equation (7)) induced by this Christoffel symbols has only one independent component,

$$R_{212}^1=R_{121}^2=\frac{{\left({\mathsf{\partial }}_{2}n\right)}^2}{n^2}+\frac{{\left({\mathsf{\partial }}_{1}n\right)}^2}{n^2}-\frac{{\mathsf{\partial }}_2^{2}n}{n}-\frac{{\mathsf{\partial }}_1^{2}n}{n}$$

(14)

inasmuch as $R_{111}^1=R_{112}^1=R_{121}^1=R_{122}^1=R_{211}^1=R_{222}^1=R_{111}^2=R_{122}^2=R_{211}^2=R_{212}^2=R_{221}^2=R_{222}^2=0$ and $R_{221}^1=R_{112}^2=-R_{212}^1$ by skew-symmetry.

Now, due to the fact that the refraction index n is a smooth function (Definition 4), $n^{-1}$ is so, so it induces the metric tensor (Lemma 1)

$${\mathbf{g}}^{\prime }=n^{-2}\delta =\frac{1}{n^2}\delta ={\mathbf{g}}^{-1}$$

where the last equality follows because ${\mathbf{g}}^{\prime }$ is the inverse of the metric tensor $\mathbf{g}$ of seismic space $\left(\mathbb{M},\mathbf{g};n\right)$. ${\mathbf{g}}^{\prime }$ induces the following Christoffel symbols,

$${{\Gamma }^{\prime }}_{11}^1=-\frac{{\mathsf{\partial }}_{1}n}{n}{{\Gamma }^{\prime }}_{12}^1={{\Gamma }^{\prime }}_{21}^1=-\frac{{\mathsf{\partial }}_{2}n}{n}{{\Gamma }^{\prime }}_{22}^1=\frac{{\mathsf{\partial }}_{1}n}{n}$$

$${{\Gamma }^{\prime }}_{11}^2=\frac{{\mathsf{\partial }}_{2}n}{n}{{\Gamma }^{\prime }}_{12}^2={{\Gamma }^{\prime }}_{21}^2=-\frac{{\mathsf{\partial }}_{1}n}{n}{{\Gamma }^{\prime }}_{22}^2=-\frac{{\mathsf{\partial }}_{2}n}{n}$$

and consequently, the following is the only independent component of curvature tensor (7)

$${R^{\prime }}_{212}^1={R^{\prime }}_{121}^2=-\frac{{\left({\mathsf{\partial }}_{2}n\right)}^2}{n^2}+\frac{{\mathsf{\partial }}_2^{2}n}{n}-\frac{{\left({\mathsf{\partial }}_{1}n\right)}^2}{n^2}+\frac{{\mathsf{\partial }}_1^{2}n}{n}$$

(15)

since by symmetry, ${R^{\prime }}_{111}^1={R^{\prime }}_{112}^1={R^{\prime }}_{121}^1={R^{\prime }}_{122}^1={R^{\prime }}_{211}^1={R^{\prime }}_{222}^1={R^{\prime }}_{111}^2={R^{\prime }}_{122}^2={R^{\prime }}_{211}^2={R^{\prime }}_{212}^2={R^{\prime }}_{221}^2={R^{\prime }}_{222}^2=0$ and by anti-symmetry, ${R^{\prime }}_{221}^1={R^{\prime }}_{112}^2=-{R^{\prime }}_{212}^1$.

By inspection to Equations (14) and (15), it can be observed that Equation (13) holds $\forall k,l,m,n=1,2$. □

Corollary 2.

By Theorem 2, a geodesic ${\kappa }^{\prime }$ (satisfying Equation (12)) in the reverse seismic space $\left(\mathbb{M},{\mathbf{g}}^{-1};n^{-1}\right)$ is a seismic ray.

For a given seismic space, $\left(\mathbb{M},\mathbf{g};n\right)$, the seismic rays ${\kappa }^{\prime }$ in their corresponding inverse seismic space $\left(\mathbb{M},{\mathbf{g}}^{-1};n^{-1}\right)$ follow the laws of reflection.

3.2. Seismic Ray Tracing Algorithmics

In summary, the steps to determine the trajectories of seismic rays in this formalism are:

1.

Determine the velocity model $v=v\left(q^i\right)$ of the region under study, and with it, the refractive index n associated (see Equation (10)).

2.

Construct the seismic space $\left({\mathbb{M}}^d,\mathbf{g};n\right)$, by constructing its metric tensor $\mathbf{g}$ induced by n.

3.

Obtain the Christoffel symbols associated with said metric tensor.

4.

Determine the geodesic (differential) equations that would correspond to the seismic rays in said means (Equations (6) and (12)).

5.

Establish the conditions of the seismic source (initial conditions for the geodesic equations).

6.

Resolution of the geodesic equations.

Regarding the last point, the geodesic equations (Equation (6)) yield the parametric geodesic trajectories (i.e., $q^i=q^i\left(\lambda \right)$), and these are usually a set of (if $i=1,2,3$) coupled second-order differential equations, which usually do not have analytic solution for cases where the velocity model is realistic. Thus, they are usually solved via numerical methods.

4. Results (Geophysical Applications)

For the purpose of exemplifying the usage of the theory herein developed in specific geophysical scenarios, in this section simplified cases are presented, for which 2D examples are chosen to simplify its visualisation, i.e., those in which $v=v(x,z)$.

4.1. Two Different Homogeneous Strata Separated at z = 2500 km

The strata interface with the following analytic speed model is simulated,

$$v(x,z)=10+3arctan\left(2500-z\right),$$

(16)

which yields the following geodesic equations,

$$\ddot{x}\left(t\right)=-\frac{6\dot{x}\left(t\right)\dot{z}\left(t\right)}{\left({(2500-z\left(t\right))}^2+1\right)\left(3{tan}^{-1}(2500-z\left(t\right))+10\right)},$$

(17)

and

$$\ddot{z}\left(t\right)=\frac{3\left(\dot{x}{\left(t\right)}^2-\dot{z}{\left(t\right)}^2\right)}{\left({(2500-z\left(t\right))}^2+1\right)\left(3{tan}^{-1}(2500-z\left(t\right))+10\right)},$$

(18)

which yields the ray diagram in the geophysical strata shown in Figure 1, with initial conditions $x\left(0\right)=z\left(0\right)=0$ km, $v_x\left(0\right)=6$ km/s and $v_z\left(0\right)=5$ km/s.

4.2. Simplified Speed Profile of Earth’s Mantle

In order to try to emulate the behaviour of P-wave speed along the mantle, the velocity model is simplified by neglecting the well-known 410 and 660 discontinuities, and before reaching the D″ layer [58]. The simplified velocity model for the mantle herein applied is given by

$$v(x,z)=6{\left(1+z\right)}^{1/9}$$

(19)

which yields the following geodesic (refraction) equations:

$$\ddot{x}\left(t\right)=\frac{2\dot{x}\left(t\right)\dot{z}\left(t\right)}{9\left[1+z\left(t\right)\right]},$$

(20)

and

$$\ddot{z}\left(t\right)=\frac{{\dot{z}}^2\left(t\right)-{\dot{x}}^2\left(t\right)}{9\left[1+z\left(t\right)\right]}.$$

(21)

The obtained ray diagram in the geophysical strata shown in Figure 2. The three rays $R_1,R_2$ and $R_3$ depicted in such a figure have the same initial position conditions, $x\left(0\right)=0$ and $z\left(0\right)=0$, as well as $v_x\left(0\right)=50$ km/s and $v_z\left(0\right)=283.56,137.37,86.60$ km/s, respectively, as speed initial conditions.

In this case, as the velocity model in Equation (19) is simple enough to obtain the Ricci scalar,

$$R=R\left(z\right)=-\frac{8}{{(1+z)}^{16/9}}$$

(22)

which indicates that for the underground, the curvature is negative.

4.3. Modelling of Irregular Interfaces between Strata

It is indeed well known that real life features between different geological bodies are not as regular as those shown in Section 4.1 or as smooth as those shown in Section 4.2. For a correct modelling of these irregularities, an interpolation of a function to delineate realistically the interfaces could yield the best results. Nevertheless, in order to show that irregularities change dramatically the trajectories of reflected and refracted rays, in this section the model presented in Section 4.1 is modified with a sawtooth function.

Suppose an irregular horizontal interface between two strata with different constant P-wave speed. In order to model it and thus, to show the effect of adjacent initial conditions geodesics in an interface with different flanks, the interface is modeled as a sawtooth function. The corresponding velocity model would be given by:

$$v(x,y)=10-3arctan\left(arctan\left(2\left(\frac{1}{4}sin\left(\frac{x\left(t\right)}{2}\right)+\frac{1}{2}sin\left(\frac{x\left(t\right)}{4}\right)-sin\left(\frac{x\left(t\right)}{8}\right)-\frac{1}{3}sin\left(\frac{3x\left(t\right)}{8}\right)-\frac{1}{5}sin\left(\frac{5x\left(t\right)}{8}\right)\right)+z\left(t\right)+100\right)\right).$$

(23)

In Figure 3, two different ray trajectories for slightly different initial conditions are shown. Both rays in both figures come out with $v_x=10$ m/s and $v_z=10tan(50°),10tan(51°)$ m/s. In Figure 3a: $x_0=400$ km and $z_0=0$ km. In Figure 3b: $x_0=370$ km and $z_0=0$ km.

4.4. Geodesics through a Salt Dome

It is well known in seismics and in geophysical exploration that the presence of salt domes yield to blurriness in the imaging of subsalt regions. This is due to the high speed contrast between the salt body and its surroundings. Thus, a velocity model for a realistic geophysical model which considers a salt dome located at Perdido’s Fold Belt is considered here, and is given by:

$$v(x,z)=\left(1000-\frac{v_1(x,z)}{{10}^{13}}\right)+v_2(x,z);$$

(24)

where

$$v_1(x,z)=x\left(t\right)p_{01}+x{\left(t\right)}^2{p}_{02}+x{\left(t\right)}^3{p}_{03}+x{\left(t\right)}^4{p}_{04}+x{\left(t\right)}^5{p}_{05}+z\left(t\right)p_{10}+$$

$$x\left(t\right)z\left(t\right)p_{11}+x{\left(t\right)}^{2}z\left(t\right)p_{12}+x{\left(t\right)}^{3}z\left(t\right)p_{13}+x{\left(t\right)}^{4}z\left(t\right)p_{14}+z{\left(t\right)}^2{p}_{20}+$$

$$x\left(t\right)z{\left(t\right)}^2{p}_{21}+x{\left(t\right)}^{2}z{\left(t\right)}^2{p}_{22}+x{\left(t\right)}^{3}z{\left(t\right)}^2{p}_{23}+z{\left(t\right)}^3{p}_{30}+x\left(t\right)z{\left(t\right)}^3{p}_{31}+$$

$$x{\left(t\right)}^{2}z{\left(t\right)}^3{p}_{32}+z{\left(t\right)}^4{p}_{40}+x\left(t\right)z{\left(t\right)}^4{p}_{41}+z{\left(t\right)}^5{p}_{50}+p_{00},$$

(25)

and

$$v_2(x,z)=4000exp\left(-\left(\frac{{\left(z\left(t\right)-x_0\right)}^4}{2a_x^2}+\frac{{\left(x\left(t\right)-y_0\right)}^4}{2a_y^2}\right)\right),$$

(26)

where $p_{00}=2569$, $p_{10}=-165.2$, $p_{01}=-0.9672$, $p_{20}=1113$, $p_{11}=59.17$, $p_{02}=$ $-15.96$, $p_{30}=1630$, $p_{21}=-53.88$, $p_{12}=34.73$, $p_{03}=36.83$, $p_{40}=-373.5$, $p_{31}=-8.61$, $p_{22}=-5.29$, $p_{13}=-14.09$, $p_{04}=5.541$, $p_{50}=-470.2$, $p_{41}=23.28$, $p_{32}=-15.3$, $p_{23}=-0.7027$, $p_{14}=-0.9127$, $p_{05}=-12.62$, $x_0=172$, $y_0=483$, $a_x=500$, and $a_y=50000$.

The latter equations set a domain of $0\le x\le 900$ km by $0\le z\le 600$ km. The velocity model yields to the geodesic equations which is not here presented due to their extension. In order to assess the behaviour of geodesics in this seismic space, three experiments were performed as follows: a source of seismic rays at ground level ($z=0$) is located, in both the left and the right corners of the domain ($x=0$ and $x=900$), as well as at the centre of the domain ($x=450$). In these locations of the source, the geodesic equations were solved several times by thoroughly sweeping the space of initial conditions, throwing geodesics each $\pi /36$ rad ($5^{\circ }$). The results of refraction and reflection experiments are shown in Figure 4 and Figure 5, respectively.

5. Discussion

5.1. Assessment of the Ray Tracing in the Previous Geophysical Setups

5.1.1. Two Different Homogeneous Strata Separated at $z=2500$ km

Although this is a very simple example, it serves to exemplify both the usage of this method as well as its accuracy. Firstly, it is important to note, from Figure 1, that within the sections of constant P-wave speed (red or blue regions), this formalism reproduce the straight lines which are expected within this kind of media.

According to the velocity model shown in Section 4.1, the diffraction indexes for both regions are $n_1=0.0679755$ and $n_2=0.189078$. On the other side, the refraction angles are ${\alpha }_1=50.19$ and ${\alpha }_2=16.02$. With such values, it is straightforward to verify that Snell’s law (law of refraction),

$$n_{1}sin\left({\alpha }_1\right)=n_{2}sin\left({\alpha }_2\right),$$

(27)

is fully satisfied.

On the other side, for the reflection process, it can be observed that the obtained reflection angle is given by $\beta =50.{19}^{\circ }$, which is consistent with the laws of reflection.

5.1.2. Simplified Speed Profile in Earth’s Mantle

Figure 2 reveals the effect of curvature in three rays travelling from the origin. It clearly affects more rays which begin travelling with a lower slope. Thus, the only way a ray could travel in real straight lines is those that travel fully vertically.

Although rays $R_1$ and $R_2$ would definitely impact the D” layer (at the base of the mantle), generating reflected rays which would be detected on Earth in a symmetric position with respect to the emission point (as the model is symmetric on x, i.e., only depends on the depth), it is clear that ray $R_3$ (which most departs from its original straight trajectory) deviates so much due to the sudden increase in P-wave speed in the first tens of metres that it would never reach D” layer and would be received at ground level in a nearer position from the origin than, for instance, ray $R_2$.

When equipped with simple enough synthetic models as this simplified one for the mantle, simple expressions as the Ricci scalar in Equation (22) can be obtained. As this scalar turns to be a function of coordinates, in this case $R=R\left(z\right)$, when evaluating it at different depths it would provide a number which characterises the curvature at this region. If the scalar turns to take values near zero, the curvature is low, and vice versa. Sometimes, the Ricci scalar provides a better insight than, for instance, density plots of the velocity model (as shown in Figure 2) as it shows directly the curvature undertaken by a ray.

Furthermore, even in cases where the Ricci scalar yields to exceptionally long analytic expressions, it can always be evaluated to numerical results in a computer, or better yet, plotted in the domain so as to achieve a better grasp of its behaviour. For instance, a graphic for the Ricci scalar in Equation (22) can be observed in Figure 6. As it can be observed, the curvature is much greater in the first tens of kilometres in the underground and then diminishes drastically after, let us say, 50 to 100 km. Below such depths, although not zero, the curvature is so small that seismic rays travel in almost straight lines. This effect is clearly reflected in rays that arrive at the most deep part of the section shown in Figure 2 (black line).

As for cases in which the Ricci scalar takes a much complicated expression, or for those in which $R=R(x,z)$ for instance, they can be plotted in a surface (3D Plot) graph to observe the behaviour of curvature in the full domain.

5.1.3. Modelling of Irregular Interfaces between Strata

Although the modelling of the interface between two homogeneous strata is ideal, the essence of the irregularities is fully captured by this example. As can be observed in Figure 3, a very slight difference in speed initial conditions of just one degree between adjacent emitted rays results in large differences in the reflected rays, which would be received in positions separated several kilometres at the surface. This is due to the nature of the irregularities introduced in the model at the separation of the strata, which presents different flanks in order to simulate real interfaces which are highly irregular.

This example shows the great flexibility of this tool to model highly complex interfaces between geological structures. It probably remains the most difficult task to make realistic models of analytic speed models of underground in order to find their corresponding seismic ray pattern.

5.1.4. Geodesics through a Salt Dome

It is important to remark that the model in Equation (24) is not symmetric as it intensely depends on both variables, $(x,z)$. Although this is not observable in the density plots of Figure 4 and Figure 5, the effects of asymmetry are evident in the behaviour of seismic rays in such figures. It should noted that the results of refraction and reflection are shown separately (in Figure 4 and Figure 5, respectively) for the sake of clarity.

The experiments shown in Figure 4 demonstrate the well-known effect that the high-speed contrast between the salt dome and the surroundings impedes seismic rays from penetrating the salt dome, hindering surveys below the salt body. Consequently with the results of Section 4.2, in Figure 4b where the source of seismic rays is located above the salt body, only almost vertical (refracted) rays penetrate the salt dome, crossing it and continuing to travel through the subsoil.

If the (refraction) experiment is performed by locating the source in the right corner of the seismic section (see Figure 4c), as the rays impinge with large angles (with respect to the vertical), they are all unable to penetrate the salt dome, being instead reflected and thus received at ground level.

On the other side, if the experiment is performed by locating the source in the other (left) corner, the asymmetry of the synthetic model arises. In this setup (Figure 4a), while the incident rays at large angles (with respect to the vertical) typically reflect at the top of the salt dome, there are (more vertical) rays that travel below the dome and are curved by the higher P-wave speed at the low part of the model, returning to its upper part. These latter rays impinge with the base of the salt body, some being reflected and some (the most vertical ones) refracted, crossing the dome and finally arriving at ground level with higher arrival times but closer to the source, complicating the interpretation of the seismic section.

As for the portion of the energy of each ray which is reflected (not refracted), the overall results for the same three experiments can be observed in Figure 5. The experiment where the source was located in the middle part above the salt dome is shown in Figure 5b; it can be observed that in this case, all the rays that impinge in the top part of the salt body cross it effectively, and continue travelling far below the salt dome in a notably vertical fashion, which does not return back to ground level.

As for both experiments with the source located at the left and right corners of the section (Figure 5a and Figure 5c, respectively), it can be observed that the fraction of energy of the incident rays enters the salt dome and stays enclosed within it in several internal reflections until finally (and probably very weakly) exiting the dome at the opposite part of the body. Nevertheless, in the case of Figure 5c and due to the asymmetry of the synthetic model, some rays emerge from the left side of the salt dome and are directed to the surface, which would be confused with rays coming from the refraction experiment, hindering the interpretation of the results.

As an overall conclusion for the experiments performed in Section 4.4, neither refracted nor reflected rays in any of the performed experiments are able to provide information about the base of the salt dome nor below which can be retrievable at the surface level. Thus, this formalism is consistent with traditional seismic surveys which always find it difficult to illuminate these regions. In this sense, Figure 4 and Figure 5 provide a more visual (physically and geometrically) insight into the outcome of seismic rays in the presence of bodies with high P-wave speed contrasts, like salt domes.

5.2. Strengths and Weaknesses of the Proposed Methodology

The strengths and weaknesses of using differential geometry for mapping salt domes in geophysical exploration can be condensed as follows:

5.2.1. Strengths

1.

Accurate modelling of Complex Geometries.

• Precise Representation: Differential geometry allows for the accurate modelling of the intricate shapes and structures of salt domes, which are often highly irregular and complex. This precision is crucial for understanding the subsurface environment.
• Curved Space Analysis: The use of differential geometry facilitates the analysis of curved spaces, providing a more realistic representation of geological formations compared to traditional Euclidean methods.

2.

Enhanced Data Interpretation.

• Detailed Structural Insights: By employing differential geometry, geoscientists can gain deeper insights into the structural complexities of salt domes, improving the interpretation of seismic data.
• Geophysical Data Integration: Differential geometry provides a robust framework for integrating various types of geophysical data, such as seismic, gravity, and magnetic data, leading to a more comprehensive understanding of the subsurface.

3.

Improved Seismic Imaging.

• Advanced Algorithms: Differential geometry enables the development of advanced algorithms for seismic imaging, such as those used in ray tracing and wave propagation modelling. These algorithms can more accurately predict the paths of seismic waves in complex media.
• Enhanced Resolution: The application of differential geometric methods can lead to enhanced resolution in seismic imaging, allowing for better identification and delineation of salt domes and associated structures.

4.

Theoretical Foundations.

• Rigorous Mathematical Framework: Differential geometry provides a rigorous mathematical framework for modelling and analysing geological structures. This foundation ensures that the methods used are theoretically sound and can be systematically improved and validated.

5.2.2. Weaknesses

1.

Computational Complexity.

• High Computational Demand: The application of differential geometry to seismic data processing and interpretation is computationally intensive. The algorithms involved often require significant computational resources and time, which can be a limitation in practical scenarios.
• Algorithm Complexity: Developing and implementing differential geometric algorithms can be complex, requiring specialised knowledge and expertise in both mathematics and geophysics.

2.

Data Quality and Availability.

• Dependence on High-Quality Data: The effectiveness of differential geometric methods is highly dependent on the quality and resolution of the input data. Poor data quality can lead to inaccurate models and interpretations.
• Data Integration Challenges: Integrating different types of geophysical data using differential geometry can be challenging, especially when dealing with disparate data sources and varying resolutions.

3.

Practical Implementation.

• Specialised Training: Geoscientists and engineers need specialised training to effectively use differential geometric methods. This requirement can be a barrier to widespread adoption in the industry.
• Cost Considerations: Implementing differential geometric techniques can be costly due to the need for high-performance computing resources and specialised software.

4.

Approximation and Simplification.

• Model Simplifications: Despite its precision, differential geometry often involves approximations and simplifications to make the problem tractable. These simplifications can sometimes overlook finer details of the geological structures.
• Assumptions and Limitations: The models used in differential geometry may rely on certain assumptions that do not always hold true in all geological settings, potentially leading to errors in interpretation.

5.3. Comparison between Techniques in Geophysical Exploration

In the context of geophysical exploration, particularly in mapping salt domes, both differential geometric methods and traditional Euclidean techniques like ray tracing using the Eikonal equation have their unique advantages and limitations. To provide a comprehensive understanding, it is essential to compare these approaches, highlighting their respective applications, strengths, and weaknesses, and relating them to relevant studies in the field.

5.3.1. Differential Geometry

These methods offer a robust framework for dealing with the intrinsic curvature of geological formations. Its application extends to modelling complex structures like salt domes, which are challenging to represent accurately with Euclidean methods.

Strengths

• Curved Space Analysis. Differential geometry can naturally handle the curvature and complexity of geological structures, leading to more accurate models.
• Advanced Seismic Imaging. Techniques derived from differential geometry, such as Riemannian wavefield extrapolation, improve seismic imaging by accounting for the curvature and heterogeneity of the subsurface.

Weaknesses

• Computationally Intensive. These methods often require significant computational resources, making them less practical for real-time applications.
• Complex Implementation. The mathematical complexity of differential geometry necessitates specialised expertise and training.

5.3.2. Euclidean Methods

Traditional Euclidean methods, including ray tracing based on the Eikonal equation, have been widely used in geophysical exploration due to their relative simplicity and computational efficiency.

Strengths

• Simplicity and Speed. Euclidean methods are generally simpler to implement and computationally less demanding, making them suitable for real-time applications.
• Widely Established. These methods are well-established and widely used in the industry, with extensive resources and tools available for practitioners.

Weaknesses

• Limited by Simplifications. Euclidean methods often assume straight-line propagation of seismic waves, which can lead to inaccuracies in complex geological settings.
• Less Accurate in Curved Spaces. These methods struggle with accurately modelling wave propagation in highly curved and heterogeneous environments, such as salt domes.

5.3.3. Comparative Studies

Some studies have compared the effectiveness of differential geometric methods and Euclidean techniques in geophysical imaging and subsurface modelling.

For instance, Ref. [59] finds that Riemannian wavefield extrapolation provides more accurate seismic images in complex geological settings, such as salt domes, by better accounting for the curvature of the subsurface. In contrast, Eikonal-based ray tracing, while faster, showed limitations in accurately representing wave paths in these environments.

Also, Ref. [60] uses finite difference methods on curved grids, inspired by differential geometry, and demonstrated improved accuracy in modelling wave propagation through complex structures compared to traditional Euclidean grid methods. This approach reduced numerical artefacts and better captured the effects of geological heterogeneities.

On the other side, Ref. [61] compares the geometric optics approaches in Riemannian and Euclidean spaces, concluding that the Riemannian approach offers superior resolution and accuracy in complex subsurface models. The Euclidean methods, while computationally efficient, often missed finer details in areas with significant curvature.

The choice between differential geometric methods and Euclidean techniques in geophysical exploration depends on the specific requirements of the application. In Section 5.3.1, it is explained that differential geometry provides enhanced accuracy and detailed modelling capabilities, making it suitable for complex geological structures like salt domes. However, it comes with higher computational costs and implementation complexity. On the other hand, Euclidean methods (see Section 5.3.2) offer simplicity and speed, beneficial for real-time applications but may fall short in accurately representing complex subsurface features.

Integrating insights from both approaches, leveraging the strengths of differential geometry for detailed modelling and Euclidean methods for computational efficiency, can provide a balanced and effective strategy in geophysical exploration. Continued comparative studies and advancements in computational techniques will further enhance the capabilities and applications of these methods in the field.

6. Conclusions

This work presents the practical implementation of differential geometry for calculating seismic ray trajectories in non-homogeneous sections. The core innovations include:

1.

Geodesic Trajectory Calculation: Using differential geometry to model seismic ray paths as geodesics in a pseudo-Riemannian metric, derived from P-wave velocities.

2.

Application to Complex Geological Setups: Demonstrating the method’s effectiveness in various geological configurations, including idealised models and salt dome structures.

3.

Refraction and Reflection Analysis: Providing a clear methodology for understanding seismic wave behaviour in complex media.

While promising, this method has limitations. It requires an analytic velocity model, necessitating prior knowledge of the region’s geology. Additionally, the method does not estimate the amount of reflected and refracted energy or attenuation through media. Future developments could focus on:

1.

Refinement of Velocity Models: Improving the accuracy and realism of analytic models to enhance the precision of ray path predictions.

2.

Energy Attenuation and Reflection Estimations: Developing algorithms to estimate energy distribution and attenuation within the media.

3.

Integration with Seismic Data: Combining this method with traditional seismic prospecting techniques for more comprehensive geological interpretations.

Finally, the use of differential geometry in seismic ray tracing offers significant advancements in modelling and interpreting geological features. Future research should aim to address the current limitations and explore broader applications in geophysical exploration.

Future Work

The application of differential geometry in geophysical (and in particular, in oil) exploration is a rapidly evolving field with numerous opportunities for future research and development. Several key areas hold potential for significant advancements:

1.

Enhanced Computational Techniques.

(a)

Development of Efficient Algorithms:

• Adaptive Algorithms: Future research can focus on developing adaptive algorithms that adjust computational resources based on the complexity of the geological structure being modelled. These algorithms would improve efficiency without compromising accuracy.
• Parallel Computing: Leveraging parallel computing and GPU acceleration can significantly reduce the computational time required for differential geometric methods, making them more practical for real-time applications.

(b)

Machine Learning Integration:

• Hybrid Models: Combining differential geometric methods with machine learning algorithms can enhance the prediction and interpretation of complex subsurface structures. Machine learning can help identify patterns and features that are difficult to capture with traditional methods alone.
• Data Augmentation: Machine learning techniques can be used to generate synthetic data that complements existing datasets, improving the robustness and accuracy of differential geometric models.

2.

Improved Data Acquisition and Processing.

(a)

High-Resolution Data Collection:

• Advanced Sensors: The development of advanced sensors and data acquisition techniques, such as 4D seismic monitoring, can provide higher resolution data that is more suitable for differential geometric analysis.
• Integrated Geophysical Surveys: Conducting integrated geophysical surveys that combine seismic, gravity, magnetic, and electromagnetic methods can provide a comprehensive dataset for more accurate modelling of geological structures.

(b)

Data Fusion Techniques:

• Multi-source Data Integration: Research into techniques for effectively integrating data from multiple sources (e.g., seismic, well logs, and satellite imagery) can enhance the accuracy of differential geometric models.
• Uncertainty Quantification: Developing methods for quantifying and reducing uncertainties in integrated datasets will be crucial for improving the reliability of differential geometric analyses.

3.

Advanced modelling and Simulation.

(a)

Complex Geological Structures:

• Subsurface Heterogeneity: Future work can focus on improving the ability of differential geometric methods to model highly heterogeneous subsurface environments, such as those with significant variations in rock properties.
• Dynamic modelling: Developing dynamic models that account for changes in geological structures over time, such as those caused by tectonic activity or fluid extraction, can provide more accurate predictions for petroleum exploration.

(b)

Incorporating Physical Processes:

• Wave-Propagation modelling: Further research can enhance the modelling of wave propagation in complex media by incorporating additional physical processes, such as anisotropy and viscoelasticity, into differential geometric frameworks.
• Thermal and Mechanical Properties: Integrating thermal and mechanical property variations into differential geometric models can improve the understanding of how these factors influence seismic wave behaviour and reservoir dynamics.

4.

Field Applications and Validation.

(a)

Case Studies and Field Trials:

• Real-World Applications: Conducting case studies and field trials to validate differential geometric methods in various geological settings, such as offshore oil fields and complex onshore basins, will be essential for demonstrating their practical value.
• Benchmarking Studies: Comparing the performance of differential geometric methods with traditional techniques in real-world scenarios can help identify best practices and areas for improvement.

(b)

Industry Collaboration:

• Collaborative Projects: Encouraging collaboration between academia, industry, and government agencies can facilitate the development and implementation of differential geometric methods in petroleum exploration.
• Standardisation: Establishing industry standards for the application of differential geometry in geophysical modelling can promote wider adoption and ensure consistency in practice.