Principe du moteur

Cette page décrit les équations résolues par le moteur de simulation d'IsoFind, le schéma numérique utilisé et les conditions dans lesquelles les résultats sont fiables. Elle s'adresse aux utilisateurs qui veulent comprendre ce qui se passe sous le capot pour mieux interpréter les sorties et reconnaître les situations où le moteur atteint ses limites.

Équation résolue

Le cœur du moteur est un modèle ADE (Advection-Dispersion-Reaction) 3D simplifié qui résout, en chaque point de la grille, l'équation suivante :

            ∂C/∂t = −v · ∇C + D · ∇²C − R(C, conditions)
        

Avec C la concentration, v le champ de vitesse effective (Darcy / porosité), D le tenseur de dispersion, et R le terme réactionnel qui dépend des conditions locales interpolées depuis les échantillons. Le moteur distingue explicitement deux modes selon la nature du contaminant simulé.

Mode	Contaminants	Terme réactionnel R
Élément	Cr, Fe, Pb, Sb, As, U, Se...	Spéciation ML, fractionnement de Rayleigh, adsorption sur solide
Molécule	Pesticides, PFAS, HAP, solvants chlorés, médicaments	Dégradation premier ordre + sorption linéaire Kd = Koc × foc

Ces deux modes partagent le même schéma numérique mais mobilisent des bridges Nexus distincts : le bridge géochimique pour les éléments (spéciation redox, fractionnement), le bridge CSIA pour les molécules (voie de dégradation dominante, epsilon, demi-vie). Les chaînes parent-métabolite sont gérées dans le mode molécule par cascade automatique lorsque la base moléculaire déclare les métabolites.

Discrétisation spatiale et temporelle

Le moteur utilise un schéma aux différences finies sur la grille régulière définie dans la configuration de la scène. Les opérateurs d'advection et de dispersion sont discrétisés par des différences centrées au second ordre, avec un traitement upwind pour les cas d'advection dominante afin de limiter les oscillations numériques.

La marche en temps est explicite par défaut, ce qui impose une condition de stabilité liant le pas de temps et la taille de maille (nombres de Courant et de Péclet). Un mode implicite est disponible pour les simulations très lentes où le pas explicite devient prohibitif.

Critère	Expression	Seuil de stabilité
Courant (advection)	Co = v · Δt / Δx	< 1 pour schéma explicite
Péclet (dispersion vs advection)	Pe = v · Δx / D	< 2 pour précision
Courbure de la concentration	Neumann D · Δt / Δx²	< 0,5 pour schéma explicite

Le moteur calcule automatiquement ces critères avant de lancer la simulation et propose un pas de temps compatible. L'utilisateur peut passer outre mais reçoit un avertissement si les critères de stabilité ne sont pas respectés.

Adsorption et retard

L'adsorption est modélisée par un coefficient de retard R appliqué à la vitesse effective : la vitesse de migration du composé devient v/R au lieu de v. Le facteur R est calculé à partir du Kd local, de la densité sèche et de la porosité.

            R = 1 + (ρ_sec / θ) · Kd
        

Pour les composés organiques, Kd est dérivé du Koc tabulé dans la base de référence et de la teneur en matière organique foc locale : Kd = Koc · foc. Cette dérivation est automatique si la base de référence contient Koc pour le composé et si les couches ont une teneur en matière organique renseignée.

Dégradation et chaînes parent-métabolite

La dégradation suit une cinétique de premier ordre par défaut : dC/dt = −k · C. Lorsqu'une chaîne parent-métabolite est déclarée, la concentration du métabolite augmente du taux de dégradation du parent, pondéré par un coefficient stœchiométrique. Les chaînes cascadées (parent → métabolite → métabolite secondaire) sont résolues simultanément.

La constante cinétique k peut être uniforme dans tout le volume, ou dépendre localement des conditions. Dans le mode couplé aux conditions géochimiques, k est calculé à chaque maille à partir de la voie dominante et des conditions pH, Eh, température, oxygène dissous interpolées sur la maille.

Calcul de la signature isotopique

Le moteur propage la signature isotopique du composé en parallèle de la concentration, en utilisant l'équation de Rayleigh. À chaque pas de temps et dans chaque maille où la dégradation est active, la signature résiduelle évolue selon le fractionnement associé à la voie dominante.

            δ(t+Δt) = δ(t) + ε · ln(C(t+Δt) / C(t))
        

Cette propagation conjointe concentration + signature est ce qui rend la simulation particulièrement utile dans IsoFind par rapport à un simple modèle de transport classique : le résultat final contient une carte spatiale des signatures isotopiques attendues, directement comparable aux mesures CSIA du terrain.

Conditions aux limites

Les bords du domaine de simulation peuvent être configurés selon trois types de conditions.

Type	Description	Usage typique
Dirichlet	Concentration imposée	Entrée d'une nappe amont connue, bord de plan d'eau
Neumann	Flux imposé (parfois zéro)	Frontière étanche, bord de bassin versant
Cauchy (mixte)	Combinaison flux-concentration	Frontière perméable avec gradient

Les conditions par défaut sont Neumann nul sur les bords latéraux et le fond, Cauchy à la sortie aval, Dirichlet au niveau des sources ponctuelles. Ce jeu couvre la majorité des cas d'hydrogéologie classique.

Performance et parallélisation

Le calcul tire parti des instructions vectorisées SIMD du processeur pour les opérations spatiales, et d'un parallélisme multithread pour les tranches de la grille. La performance typique sur une grille 100 × 100 × 40 est de l'ordre de quelques centaines de pas de temps par seconde sur un poste de travail récent, soit une simulation complète de plusieurs années en moins d'une minute.

Pour les grilles très fines ou les simulations de longue durée, un mode calcul différé permet de lancer la simulation en arrière-plan pendant que l'utilisateur continue à travailler dans l'interface. Les résultats sont chargés automatiquement dès qu'ils sont disponibles.

Limites du schéma numérique

Le schéma explicite introduit une diffusion numérique artificielle lorsque le nombre de Courant approche 1. Cette diffusion peut masquer un gradient réel fort en bordure de source.
L'upwind d'ordre 1 utilisé en zone d'advection dominante lisse les fronts de concentration. Pour les cas où un front raide est critique, le passage à un schéma d'ordre supérieur est possible mais plus coûteux.
Le couplage concentration + isotopie suppose une voie dominante unique. En cas de plusieurs voies concurrentes avec des epsilons très différents, la signature simulée est une moyenne pondérée qui peut mal représenter l'effet réel.
Les milieux fortement anisotropes ou à fractures discrètes demandent une représentation explicite des hétérogénéités, que le schéma maille régulière simule difficilement.

Pour aller plus loin

Paramètres hydrogéologiques : comment les conditions d'écoulement sont renseignées.
Spéciation et propagation : couplage réactif et cinétique.
Fractionnements isotopiques : le moteur Nexus qui fournit les epsilons.