En exploration géophysique, la magnétotellurique (MT) est une technique clé. Le calcul de modélisation directe pour les modes TE (transverse électrique) et TM (transverse magnétique) en constitue un aspect fondamental.

Principes des modes TE et TM

La magnétotellurique exploite des champs électrométalliques naturels et alternatifs pour sonder la structure électrique souterraine. Comme le milieu souterrain présente généralement une hétérogénéité bidimensionnelle ou tridimensionnelle, il est nécessaire de calculer séparément les réponses des deux modes de polarisation TE et TM pour une compréhension complète de la structure géologique.

En pratique, en raison de l'hétérogénéité du milieu souterrain, les courbes de résistivité apparente pour les modes TE et TM présentent souvent des différences significatives. Comprendre les caractéristiques de réponse de ces deux modes est essentiel pour une interprétation correcte des structures géologiques réelles.

Fondements mathématiques de la modélisation directe

Équations de contrôle

En partant des équations de Maxwell, on peut dériver le problème aux limites des équations aux dérivées partielles gouvernant la modélisation directe magnétotellurique 2D.

Pour le mode TE, l'équation de contrôle s'écrit :

\( \nabla \cdot (\mu^{-1} \nabla E_y) + i\omega \sigma E_y = 0 \)

Pour le mode TM, l'équation est :

\( \nabla \cdot (\sigma^{-1} \nabla H_y) + i\omega \mu H_y = 0 \)

Où \( E_y \) et \( H_y \) sont respectivement les composantes du champ électrique et magnétique perpendiculaires au profil, \( \sigma \) est la conductivité électrique, \( \mu \) est la perméabilité magnétique et \( \omega \) est la pulsation.

Méthodes de discrétisation numérique

La modélisatino directe magnétotellurique repose principalement sur des méthodes de simulation numérique. Les approches courantes sont :

• Méthode des éléments finis (MEF) : adaptée aux terrains complexes et aux conditions géologiques variées.
• Méthode des différences finies (MDF) : relativement simple à implémenter.
• Méthode des éléments spectraux (SEM) : offre une précision élevée et une convergence rapide.

Parmi celles-ci, la méthode des éléments finis est la plus répandue. Pour améliorer la précision, on peut utiliser des techniques d'interpolation bilinéaire ou bi-quadratique.

Implémentation en Python

Ci-dessous, un cadre d'implémentation en Python basé sur la méthode des éléments finis pour les modélisations directes des modes TE et TM.

Structure principale du programme

import numpy as np
from scipy.sparse import lil_matrix, csc_matrix
from scipy.sparse.linalg import spsolve

def magneto_forward_model(config, frequency_list):
    """
    Calcul de modélisation directe 2D pour les modes TE et TM.
    
    Paramètres:
        config (dict): Dictionnaire contenant le modèle de conductivité,
                       le maillage et les paramètres.
        frequency_list (array): Liste des fréquences à calculer.
    
    Renvoie:
        rho_te (ndarray): Résistivité apparente pour le mode TE.
        rho_tm (ndarray): Résistivité apparente pour le mode TM.
    """
    n_f = len(frequency_list)
    n_stations = len(config['station_coords'])
    rho_te = np.zeros((n_stations, n_f))
    rho_tm = np.zeros((n_stations, n_f))
    
    for idx, freq in enumerate(frequency_list):
        pulsation = 2 * np.pi * freq
        
        # Calcul pour le mode TE
        e_field_te = solve_te_mode(config, pulsation)
        
        # Calcul pour le mode TM
        h_field_tm = solve_tm_mode(config, pulsation)
        
        # Calcul des résistivités apparentes
        rho_te[:, idx] = apparent_resistivity_te(e_field_te, config, pulsation)
        rho_tm[:, idx] = apparent_resistivity_tm(h_field_tm, config, pulsation)
    
    return rho_te, rho_tm

Solveur pour le mode TE

def solve_te_mode(cfg, omega):
    """Résolution du système d'équations pour le mode TE."""
    
    nodes, elems = create_mesh(cfg['domain'], cfg['mesh_settings'])
    n_nodes = nodes.shape[0]
    
    # Initialisation des matrices creuses
    stiffness = lil_matrix((n_nodes, n_nodes))
    mass = lil_matrix((n_nodes, n_nodes))
    
    for elem in elems:
        node_indices = elem
        local_coords = nodes[node_indices]
        sigma = element_conductivity(cfg['model'], node_indices)
        
        # Calcul des matrices locales
        K_local, M_local = te_element_matrices(local_coords, sigma, omega)
        
        # Assemblage
        for i_local, i_global in enumerate(node_indices):
            for j_local, j_global in enumerate(node_indices):
                stiffness[i_global, j_global] += K_local[i_local, j_local]
                mass[i_global, j_global] += M_local[i_local, j_local]
    
    # Application des conditions limites
    stiffness, mass = apply_te_bc(stiffness, mass, nodes, cfg)
    
    # Second membre (source)
    rhs = te_source_vector(nodes, cfg)
    
    # Résolution du système linéaire: (K + i*ω*μ*M) * E = source
    system_matrix = stiffness + 1j * omega * cfg['mu'] * mass
    system_matrix = csc_matrix(system_matrix)
    
    solution = spsolve(system_matrix, rhs)
    return solution

Solveur pour le mode TM

def solve_tm_mode(cfg, omega):
    """Résolution du système d'équations pour le mode TM."""
    
    nodes, elems = create_mesh(cfg['domain'], cfg['mesh_settings'])
    n_nodes = nodes.shape[0]
    
    stiffness_tm = lil_matrix((n_nodes, n_nodes))
    mass_tm = lil_matrix((n_nodes, n_nodes))
    
    for elem in elems:
        node_indices = elem
        local_coords = nodes[node_indices]
        sigma = element_conductivity(cfg['model'], node_indices)
        rho = 1.0 / sigma  # Résistivité
        
        # Calcul des matrices locales pour le mode TM
        K_local, M_local = tm_element_matrices(local_coords, rho, omega)
        
        # Assemblage
        for i_local, i_global in enumerate(node_indices):
            for j_local, j_global in enumerate(node_indices):
                stiffness_tm[i_global, j_global] += K_local[i_local, j_local]
                mass_tm[i_global, j_global] += M_local[i_local, j_local]
    
    # Conditions limites pour le mode TM
    stiffness_tm, mass_tm = apply_tm_bc(stiffness_tm, mass_tm, nodes, cfg)
    
    # Second membre
    rhs_tm = tm_source_vector(nodes, cfg)
    
    # Résolution
    A_tm = stiffness_tm + 1j * omega * cfg['mu'] * mass_tm
    A_tm = csc_matrix(A_tm)
    
    solution = spsolve(A_tm, rhs_tm)
    return solution

Calcul des résistivités apparentes

def apparent_resistivity_te(field_e, cfg, omega):
    """Calcule la résistivité apparente pour le mode TE."""
    
    # Calcul numérique du gradient du champ E
    grad_ex = numerical_gradient(field_e, cfg['station_coords'])
    
    # Le champ magnétique H_z est obtenu via la loi de Faraday
    H_z = (-1.0 / (1j * omega * cfg['mu'])) * grad_ex
    
    # Impédance Z_xy = E_y / H_z
    Z_xy = field_e / H_z
    
    # Résistivité apparente
    rho_app = (np.abs(Z_xy)**2) / (omega * cfg['mu'])
    return rho_app

def apparent_resistivity_tm(field_h, cfg, omega):
    """Calcule la résistivité apparente pour le mode TM."""
    
    # Calcul numérique du gradient du champ H
    grad_hx = numerical_gradient(field_h, cfg['station_coords'])
    
    # Le champ électrique E_z est obtenu via la loi d'Ampère
    sigma_vals = get_station_conductivity(cfg['model'], cfg['station_coords'])
    E_z = (1.0 / sigma_vals) * grad_hx
    
    # Impédance Z_yx = E_z / H_y
    Z_yx = E_z / field_h
    
    # Résistivité apparente
    rho_app = (np.abs(Z_yx)**2) / (omega * cfg['mu'])
    return rho_app

Validation du modèle et analyse

Procédure de validation

Pour valider le code, on peut le tester avec un modèle géoélectrique de référence, comme un milieu semi-infini homogène ou une structure multicouche connue.

def validate_code():
    """Teste le code avec un modèle de référence."""
    
    # Modèle d'une couche sur un substratum
    resistivity_layer = 50    # Ω·m
    resistivity_substrate = 100 # Ω·m
    thickness_layer = 500     # m
    
    # Configuration
    test_config = {
        'model': {'layer_res': [resistivity_layer, resistivity_substrate],
                  'thicknesses': [thickness_layer]},
        'domain': {'x_range': (-2000, 2000), 'z_range': (0, 3000)},
        'station_coords': np.linspace(-1000, 1000, 21),
        'mu': 4 * np.pi * 1e-7,
        'mesh_settings': {'dx': 100, 'dz': 50}
    }
    
    # Fréquences de test
    freqs = np.logspace(-2, 2, 15)
    
    # Exécution de la modélisation
    rho_te, rho_tm = magneto_forward_model(test_config, freqs)
    
    # Solution analytique pour un modèle 1D (calcul via transformée de Hankel)
    rho_theoretical = analytical_1d_response(test_config['model'], freqs)
    
    # Comparaison graphique (pseudo-code)
    # plot_comparison(freqs, rho_te[10], rho_tm[10], rho_theoretical)

Application à des structures géologiques types

Le cadre peut être utilisé pour simuler des structures géologiques typiques telles que des graben, des horsts ou des failles. Les résultats de la modélisation directe révèlent comment les modes TE et TM répondent différemment à ces structures, fournissant une base essentielle pour l'interprétation des données réelles.