2.4.5 Modèles numériques 2D

• Méthode de l’intégrale de contour
• Exemple : gravité d’un polygone quelconque
  • Constantes physiques
  • Géométrie du problème
  • Paramètres du levé
  • Fonctions mathématiques
  • Calcul de la gravité
  • Visualisation des résultats
• Extension 3D de la méthode
• Références

---

Méthode de l’intégrale de contour

Soit un corps 2D de section quelconque dont l’extension dans la troisième dimension est considéré comme infini. On peut montrer que l’anomalie gravimétrique verticale $\Delta g_z$ produite par ce corps est égale à une intégrale de ligne autour de sa section [1] :

\[\Delta g = 2G\Delta\rho \oint z d\theta ,\]

où $G$ est la constante gravitationnelle et $\Delta\rho$ le contraste de densité entre le corps et son milieu encaissant.

Numériquement, on peut approximer le corps quelconque comme un polygone de $N$ côtés [2] et l’intégrale se simplifie à la sommation

\[\oint z d\theta = \sum_i^N Z_i ,\]

où $Z_i$ est l’intégrale de ligne pour le $i^\textrm{e}$ côté du polygone. Dans ce qui suit, les indices $i$ et ${i+1}$ réfèrent respectivement à un sommet du polygone et le prochain sommet dans le sens horaire, tel que montré à la Figure 1. Il est possible de représenter à peu près n’importe quel polygone, aussi complexe soit-il, en augmentant le nombre de côtés $N$ qui le décrit.

Figure 1. Polygone à $N=8$ côtés dont les sommets sont numérotés dans le sens horaire.

Avec les simplifications de [3], $Z_i$ est donné par

\[Z_i = A_i \left[ (\theta_i - \theta_{i+1}) +B_i \ln{\frac{r_{i+1}}{r_i}} \right] ,\]

avec les substitutions

\[A_i = \frac{(x_{i+1} - x_i)(x_i z_{i+1} - x_{i+1} z_i)}{(x_{i+1} - x_i)^2 + (z_{i+1} - z_i)^2} ,\] \[B_i = \frac{(z_{i+1} - z_i)}{(x_{i+1} - x_i)} ,\] \[\theta_i = \arctan{\left(\frac{z_i}{x_i}\right)} ,\]

et

\[r_i^2 = x_i^2 +z_i^2 .\]

---

Exemple : gravité d’un polygone quelconque

Dans cet exemple nous allons modéliser la réponse gravimétrique d’un polygone généré aléatoirement.

Constantes physiques

On définit d’abord les constantes physiques.

import numpy as np
import matplotlib.pyplot as plt

# Quelques constantes physiques
G = 6.67408e-11  # SI
rho = 1000  # kg/m3

La valeur de $\Delta\rho$ pour le polygone est arbitrairement fixée à 1000 kg/m$^3$ pour cet exemple.

Géométrie du problème

Un polygone aléatoire a été généré. Il est représenté par une suite de paires de coordonnées ($x_i, y_i$) de ses sommets. Vous pouvez importer directement les données à partir de mon site web pour reproduire l’expérience avec la fonction suivante :

import requests
from io import BytesIO

def np_load_from_url(url):
    address = BytesIO(requests.get(url).content)
    return np.loadtxt(address)

url = 'https://clberube.org/glq2200/assets/data/random-polygon.txt'
data = np_load_from_url(url)

Le polygone généré aléatoirement est celui de la Figure 2. Notez qu’il possède un grand nombre de sommets pour approximer son contour curviligne.

Figure 2. Polygone généré aléatoirement pour lequel nous allons simuler la réponse gravimétrique.

Paramètres du levé

Il faut maintenant définir les paramètres du levé gravimétrique, et écrire les équations dans une fonction.

# On définit les paramètres du levé
n_stations = 100
xs = np.linspace(-5000, 5000, n_stations)  # le gravimètre bouge sur x
zs = np.zeros(xs.shape)  # le gravimètre est à z = 0 partout

Fonctions mathématiques

On écrit l’équation (3) dans une fonction. Il faut aussi définir une fonction qui retourne la distance entre deux points, car la station de mesure n’est pas nécessairement à l’origine.

def dist(a, b):
    """Retourne la distance entre a et b.
    """
    return a - b

def integrand_talwani(x1, x2, z1, z2, density):
    """Suppose que l'intégration se fait en sens horaire.
    """     
    # Calcul des composantes
    denom = z2 - z1
    alpha = (x2 - x1) / denom
    beta = (x1 * z2 - x2 * z1) / denom  
    factor = beta / (1 + alpha * alpha)
    r1sq = (x1 * x1 + z1 * z1)  # r1^2 
    r2sq = (x2 * x2 + z2 * z2)  # r2^2
    term1 = 0.5 * (np.log(r2sq) - np.log(r1sq))  # pour éviter de calculer ln(r2^2/r1^2)
    term2 = np.arctan2(z2, x2) - np.arctan2(z1, x1)  # theta2 - theta1
    return factor * (term1 - alpha * term2)

J’ai fait quelques substitutions dans l’équation pour la rendre plus stable numériquement, à cause du terme en $\ln\left({r_{i+1}^2}/{r_i^2}\right)$. Vous pouvez faire l’exercice pour vous convaincre que cette notation est équivalente à l’équation donnée plus haut.

Calcul de la gravité

Il faut maintenant calculer la contribution de tous les côtés du polygone, en faire la sommation pour obtenir la gravité totale, et puis recommencer pour chaque station de mesure gravimétrique sur le profil $x$. Le calcul se fait ainsi :

# Loop sur les stations de mesures
gz = np.empty(n_stations)
for i in range(n_stations):      
    xp = dist(xs[i], x0)
    zp = dist(zs[i], z0)
    Z = integrand_talwani(xp, np.roll(xp, 1), zp, np.roll(zp, 1), rho)
    gz[i] = 2*G*rho*np.sum(Z)  # gravité en m/s^2

Visualisation des résultats

Finalement, on trace un graphique du polygone et de sa réponse gravimétrique (Figure 3).

if rho > 0:
    color = 'C3'  # rouge si (+)
else: 
    color = 'C0'  # bleu si (-)
fig, ax = plt.subplots(2, 1, figsize=(5, 8), sharex=True)
ax[1].set_aspect('equal')
ax[1].set_xlim((xs.min(), xs.max()))
ax[1].set_ylim((xs.min()*1.5, 0))
ax[0].plot(xs, gz, lw=1.5, color=color) 
ax[1].plot(x0, z0, lw=2, color=color)
ax[1].fill(x0, z0, alpha=0.25, color=color)
ax[0].set_ylabel('$\Delta g_z$ (m/s$^2$)')
ax[1].set_ylabel('$z$ (m)')
ax[1].set_xlabel('$x$ (m)')
plt.tight_layout()

Figure 3. Réponse gravimétrique d’un polygone quelconque offrant un contraste de densité avec son milieu encaissant.

Notez la légère asymétrie dans la réponse gravimétrique sur la Figure 3. C’est causé par l’asymétrie du polygone. Cool!

Question : D’après vous, quelle est la signification de la position en $x$ du maximum de l’anomalie?

Je vous encourage à répéter l’expérience avec votre propre polygone en définissant simplement de nouvelles coordonnées pour les sommets. Commencez avec quelque chose de simple, par exemple avec ce triangle :

x0 = [0, 2000, -2000]
z0 = [-1000, -3000, -2000]

---

Extension 3D de la méthode

D’après [2], cette méthode peut être utilisée pour modéliser la réponse de corps tridimensionnels. En effet, une fois qu’on a trouvé la gravité causée par un polygone, on peut intégrer selon la troisième dimension pour obtenir la réponse du corps. L’intégration se fait un peu comme si on voulait avoir le volume total d’une montagne à partir des lignes de contours d’une carte topographique. Il existe cependant des méthodes plus simples pour modéliser la gravité en 3D, comme celle de la discrétisation du milieu que nous verrons dans la prochaine section.

---

Références

1. Hubbert, M. K. (1948). A line-integral method of computing the gravimetric effects of two-dimensional masses. Geophysics, 13(2), 215-225.
2. Talwani, M., Worzel, J. L., & Landisman, M. (1959). Rapid gravity computations for two‐dimensional bodies with application to the Mendocino submarine fracture zone. Journal of geophysical research, 64(1), 49-59.
3. Won, I. J., & Bevis, M. (1987). Computing the gravitational and magnetic anomalies due to a polygon: Algorithms and Fortran subroutines. Geophysics, 52(2), 232-238.