En été sur la route, lorsque celle-ci est bien chauffée par le Soleil, on a parfois l’impression de voir des flaques d’eau au loin. La surface semble brillante, presque réfléchissante, alors qu’il n’y a pas d’eau sur la chaussée.
Ce phénomène est un mirage inférieur. Il ne vient pas d’une réflexion sur la route, mais d’une courbure progressive des rayons lumineux dans l’air situé juste au-dessus du sol.
La route chauffée par le Soleil réchauffe l’air proche de la chaussée. Cet air chaud est moins dense que l’air situé plus haut.
La lumière ne se propage pas de la même manière dans un air plus dense et dans un air moins dense. En effet, l’indice de réfraction de l’air varie avec sa densité.
Près de la route, l’indice est plus faible. Plus haut, l’air est plus froid, plus dense, et l’indice est plus grand. Il se forme donc un gradient vertical d’indice de réfraction.
Un rayon lumineux qui se propage vers le bas traverse des couches d’air où l’indice devient plus faible. La direction du rayon change alors progressivement. Le rayon devient de moins en moins descendant. S’il est assez rasant, il peut devenir horizontal avant de remonter ; sinon, il atteint la route.
Notre œil ne sait pas que le rayon a suivi une trajectoire courbe. Il le prolonge mentalement en ligne droite. Il croit alors que la lumière vient du sol, comme si elle avait été réfléchie par une flaque d’eau.
Objectif du modèle
On cherche maintenant à construire un modèle simple permettant de décrire cette courbure des rayons lumineux.
L’idée générale est la suivante :
- La route chaude impose un profil de température T(z).
- Ce profil de température amène un profil d’indice optique n(z).
- Le rayon lumineux obéit à une relation de conservation du type Snell-Descartes.
- Cette relation permet de prédire la trajectoire du rayon.
On travaillera dans un plan vertical.
- L’axe (Ox) est horizontal, parallèle à la route ; la coordonnée x mesure la distance le long de la route.
- L’axe (Oz) est vertical, orienté vers le haut ; la coordonnée z mesure la hauteur au-dessus de la route. La route est située à z=0.
- L’indice optique de l’air dépend seulement de la hauteur. On le notera n(z).
Contexte physique du mirage
Au-dessus d’une route chaude, la situation est la suivante :
- l’air très proche de la route est fortement chauffé ;
- cet air chaud est moins dense ;
- l’indice de réfraction de l’air chaud est légèrement plus faible ;
- l’air plus haut est plus froid, plus dense, donc d’indice légèrement plus grand.
n(z) augmente avec z.
Le gradient d’indice est l’ingrédient physique essentiel du modèle : il permet de décrire comment la direction d’un rayon évolue lorsqu’il se propage près de la route.
Modélisation de la température
Choix d’un profil simple
On cherche un modèle simple pour la température de l’air au-dessus de la route. On pose :
T(z)=T_\infty+\Delta T e^{-z/H}.Les paramètres sont :
- T_\infty : température de l’air loin de la route ;
- \Delta T : excès de température au voisinage immédiat de la route ;
- H : hauteur caractéristique de décroissance thermique.
À z=0, on a :
T(0)=T_\infty+\Delta T.Quand on monte, l’exponentielle diminue jusqu’à tendre vers 0, et :
T(z)\to T_\infty.D’où vient l’exponentielle ?
On note l’excès de température :
u(z)=T(z)-T_\infty.On suppose que cet excès décroît proportionnellement à lui-même quand on monte. Autrement dit, plus l’excès de température est grand, plus il diminue rapidement avec la hauteur.
Cela s’écrit :
\frac{du}{dz}=-\frac{u}{H}.La solution est :
u(z)=u(0)e^{-z/H}.Comme u(0)=\Delta T, on obtient :
T(z)-T_\infty=\Delta T e^{-z/H},donc :
T(z)=T_\infty+\Delta T e^{-z/H}.Ce n’est pas une loi exacte de météorologie. C’est un modèle phénoménologique, représentant la situation avec très peu de paramètres.
Dans une situation réelle, le profil dépendrait aussi, notamment :
- du vent ;
- de la turbulence ;
- de l’ensoleillement ;
- de la nature de la route ;
- de l’humidité ;
- de la convection ;
- de la distance le long de la route.
Mais pour comprendre le mécanisme du mirage, le modèle exponentiel est suffisant (et certains articles expliquent sa pertinence avec des arguments plus élaborés).
De la température à l’indice de réfraction
L’indice de l’air est très proche de 1. Dans des conditions ordinaires, on a n \simeq 1{,}0003..
La quantité intéressante est n-1. Pour un gaz, l’indice optique est lié à la densité. Dans un modèle simple, on écrit :
n-1 \propto \rho,où \rho est la masse volumique de l’air.
Avec l’équation des gaz parfaits, à pression à peu près constante sur quelques mètres, la masse volumique varie comme :
\rho \propto \frac{1}{T}.On obtient donc :
n(z)-1 \propto \frac{1}{T(z)}.Si n_\infty est l’indice de l’air loin de la route, à la température T_\infty, on peut écrire :
n(z)=1+(n_\infty-1)\frac{T_\infty}{T(z)}.En remplaçant T(z) par le modèle choisi :
n(z)=1+(n_\infty-1) \frac{T_\infty}{T_\infty+\Delta T e^{-z/H}}.Approximation de l’indice
Si \Delta T reste petit devant T_\infty, on peut recourir à une approximation très pratique.
On part de :
n(z)=1+(n_\infty-1) \frac{T_\infty}{T_\infty+\Delta T e^{-z/H}}.On réécrit le facteur dépendant de la température sous la forme :
\frac{T_\infty}{T_\infty+\Delta T e^{-z/H}} = \frac{1}{1+\frac{\Delta T}{T_\infty}e^{-z/H}}.Si \Delta T/T_\infty est petit, on utilise l’approximation :
\frac{1}{1+\varepsilon} \simeq 1-\varepsilonavec :
\varepsilon = \frac{\Delta T}{T_\infty}e^{-z/H}.On obtient :
n(z)\simeq n_\infty-\Delta n e^{-z/H},avec :
\Delta n=(n_\infty-1)\frac{\Delta T}{T_\infty}.Cette forme est très lisible :
n(z)=n_\infty-\Delta n e^{-z/H}.À z=0 :
n(0)=n_\infty-\Delta n.Quand z devient grand :
n(z)\to n_\infty.Quelques valeurs
Pour une route chaude en été, prenons :
T_\infty=303\ \mathrm{K},c’est-à-dire environ 30^\circ\mathrm{C},
\Delta T=20\ \mathrm{K},et
H=0{,}20\ \mathrm{m}.On prend aussi :
n_\infty=1{,}000265.Alors :
\Delta n=(n_\infty-1)\frac{\Delta T}{T_\infty}.Numériquement :
\Delta n \simeq 2{,}65\times10^{-4}\times\frac{20}{303}.Donc :
\Delta n\simeq 1{,}75\times10^{-5}.Le modèle approché devient :
n(z)\simeq 1{,}000265-1{,}75\times10^{-5}e^{-z/0{,}20}.Près de la route :
n(0)\simeq 1{,}0002475.Plus haut :
n(z)\to1{,}000265.La variation d’indice est très petite, de l’ordre de 10^{-5}. Pourtant, elle suffit à courber des rayons qui se propagent sur de grandes distances presque parallèlement à la route.
Représenter les profils T(z) et n(z)
Le code suivant trace les modèles du profil de température et du profil d’indice.
Voir le code Python
import numpy as np
import matplotlib.pyplot as plt
T_inf = 303.0 # K
Delta_T = 20.0 # K
H = 0.20 # m
n_inf = 1.000265
z = np.linspace(0, 1.5, 500)
T = T_inf + Delta_T * np.exp(-z / H)
n_exact = 1 + (n_inf - 1) * T_inf / T
Delta_n = (n_inf - 1) * Delta_T / T_inf
n_approx = n_inf - Delta_n * np.exp(-z / H)
fig, axes = plt.subplots(1, 2, figsize=(10, 4))
axes[0].plot(z, T - 273.15)
axes[0].set_xlabel("Hauteur z (m)")
axes[0].set_ylabel("Température T(z) (°C)")
axes[0].set_title("Température en fonction de la hauteur")
axes[0].grid(True)
axes[1].plot(z, n_exact, label="modèle issu de 1/T")
axes[1].plot(z, n_approx, "--", label="approximation exponentielle")
axes[1].set_xlabel("Hauteur z (m)")
axes[1].set_ylabel("Indice n(z)")
axes[1].set_title("Indice en fonction de la hauteur")
axes[1].grid(True)
axes[1].legend()
fig.tight_layout()
plt.show()
Loi de Snell-Descartes continue
Angle utilisé
On note \varphi l’angle que fait le rayon avec l’horizontale.
- Si le rayon (venant de la gauche) monte vers la droite, alors \varphi\gt 0.
- Si le rayon est horizontal, alors \varphi=0.
- Si le rayon descend vers la droite, alors \varphi\lt 0.
Dans un milieu stratifié, c’est-à-dire un milieu où l’indice dépend seulement de z, la loi de Snell-Descartes prend une forme continue :
n(z)\cos\varphi=C,où C est une constante le long du rayon.
Cette relation est la base du modèle géométrique du rayon. Elle signifie que lorsque le rayon traverse des zones où n(z) change, son angle \varphi change aussi pendant que le produit n(z)\cos\varphi reste constant.
Constante fixée par les conditions initiales
Si le rayon passe en z_0 avec un angle initial \varphi_0, alors :
C=n(z_0)\cos\varphi_0.Une fois C connue, on peut retrouver l’angle du rayon à n’importe quelle hauteur z par :
\cos\varphi=\frac{C}{n(z)}.Pour que cette expression ait un sens, il faut :
\frac{C}{n(z)}\leq 1.Si, en descendant, le rayon atteint une hauteur z_\ast telle que :
n(z_\ast)=C,alors :
\cos\varphi=1.Donc :
\varphi=0.Le rayon est horizontal. C’est le point de retournement.
Retrouver la trajectoire du rayon
Avant d’écrire l’équation de la trajectoire, on peut regarder localement un petit élément de rayon.
La pente du rayon est :
\frac{dz}{dx}=\tan\varphi.Or
\cos\varphi=\frac{C}{n(z)}et
\sin^2\varphi=1-\cos^2\varphi.Donc :
\sin^2\varphi= 1-\frac{C^2}{n(z)^2}.Ainsi :
\tan^2\varphi = \frac{\sin^2\varphi}{\cos^2\varphi} = \frac{n(z)^2}{C^2}-1.La trajectoire vérifie donc :
\frac{dz}{dx} = \pm\sqrt{\frac{n(z)^2}{C^2}-1}.Le signe dépend du sens vertical du rayon.
Pour un rayon qui descend vers la droite :
\frac{dz}{dx} = -\sqrt{\frac{n(z)^2}{C^2}-1}.Pour un rayon qui remonte vers la droite :
\frac{dz}{dx} = +\sqrt{\frac{n(z)^2}{C^2}-1}.On peut reconstruire la trajectoire à partir de n(z).
Solution formelle par intégration
La relation précédente donne déjà une solution formelle du problème.
En effet, pour un rayon descendant vers la droite :
\frac{dz}{dx} = -\sqrt{\frac{n(z)^2}{C^2}-1}.On peut donc écrire :
\frac{dx}{dz} = -\frac{C}{\sqrt{n(z)^2-C^2}}.Si l’on veut calculer la distance horizontale parcourue entre une hauteur initiale z_0 et une hauteur plus basse z, on obtient :
x(z) = \int_z^{z_0} \frac{C}{\sqrt{n(u)^2-C^2}}\,du.Cette formule est une solution complète sous forme d’intégrale.
Dans notre modèle approché :
n(z)=n_\infty-\Delta n e^{-z/H}.Posons, pour alléger les écritures :
a=n_\infty.On a alors :
n(z)=a-\Delta n e^{-z/H}.Comme :
\frac{dn}{dz} = \frac{a-n}{H},on peut utiliser n comme nouvelle variable d’intégration :
dz = \frac{H}{a-n}\,dn.L’intégrale devient :
\int \frac{C}{\sqrt{n^2-C^2}}\,dz = HC \int \frac{dn}{(a-n)\sqrt{n^2-C^2}}.Cette intégrale se calcule explicitement. On obtient :
I(n) = \frac{2}{\sqrt{a^2-C^2}} \operatorname{artanh} \left( \sqrt{ \frac{(a+C)(n-C)} {(a-C)(n+C)} } \right).La branche descendante, en partant de z_0, s’écrit donc :
x(z) = HC\left[I(n(z_0))-I(n(z))\right].Le point de retournement est atteint lorsque n(z_\ast)=C.
Avec le profil exponentiel, cela donne :
z_\ast = -H\ln\left(\frac{n_\infty-C}{\Delta n}\right).Au point de retournement, on a n(z_\ast)=C, donc I(n(z_\ast))=0.
La distance horizontale jusqu’au retournement vaut donc :
x_\ast = HC\,I(n(z_0)).Pour la branche montante, on repart du point de retournement et on prend l’autre signe de la racine. On obtient :
x(z) = x_\ast + HC\,I(n(z)).Cette solution est formelle et rigoureuse pour le profil approché n(z)=n_\infty-\Delta n e^{-z/H}. Elle ne concerne donc pas directement le modèle complet issu de 1/T(z), ni une atmosphère réelle plus complexe.
En pratique, elle permet de calculer x en fonction de z. Pour tracer la courbe, cela suffit : on choisit des valeurs de z, puis on calcule les valeurs correspondantes de x.
Les faibles angles seulement
La variation d’indice au-dessus d’une route chaude est très faible, typiquement de l’ordre de 10^{-5}.
Cela signifie que seuls les rayons très proches de l’horizontale sont suffisamment sensibles à cette variation pour être retournés.
Un rayon arrivant avec un angle de 45^\circ par rapport à l’horizontale sera seulement très légèrement dévié. Il ne donnera pas le mirage typique.
Pour obtenir un retournement, il faut un rayon rasant, avec un angle de l’ordre du dixième de degré ou de quelques dixièmes de degré.
C’est pourquoi les mirages de route sont observés à grande distance, avec des directions de visée presque horizontales.
Comparer trois méthodes de tracé
On va maintenant comparer trois façons de tracer les rayons, de la plus directe numériquement à la plus spécialisée pour notre profil d’indice.
La première méthode est une intégration pas à pas de type Euler. On avance le long du rayon par petits éléments ds, en mettant à jour x, z et \varphi :
\frac{dx}{ds}=\cos\varphi, \frac{dz}{ds}=\sin\varphi,et :
\frac{d\varphi}{ds} = \frac{1}{n(z)} \frac{dn}{dz} \cos\varphi.Cette méthode est volontairement élémentaire : elle est générale, mais elle dépend du pas ds et accumule une erreur numérique.
La deuxième méthode utilise la constante de Snell-Descartes pour écrire directement :
\frac{dz}{dx} = \pm\sqrt{\frac{n(z)^2}{C^2}-1}.La branche descendante utilise le signe négatif. La branche montante utilise le signe positif. On obtient alors une quadrature numérique : on connaît l’intégrale à calculer, mais on l’évalue par petits trapèzes.
Près du point de retournement, l’intégrande en z devient très grande, même si l’intégrale reste finie. Pour éviter un artefact numérique à cet endroit, on pose z=z_\ast+u^2 et on applique les trapèzes sur la variable u. Le minimum du rayon est alors traité comme un point régulier.
La troisième méthode utilise tout bonnement la formule analytique obtenue plus haut pour l’intégrale. Elle est exacte pour le profil d’indice approché :
n(z)=n_\infty-\Delta n e^{-z/H}.Elle ne prétend donc pas décrire toute l’atmosphère réelle exactement : elle donne la solution de référence pour le modèle simplifié utilisé ici.
Les trois graphes ci-dessous utilisent les mêmes paramètres, les mêmes angles initiaux et les mêmes échelles sur les axes.
Paramètres utilisés :
- T_\infty=303\ \mathrm{K} ;
- \Delta T=20\ \mathrm{K} ;
- H=0{,}20\ \mathrm{m} ;
- n_\infty=1{,}000265 ;
- z_0=1{,}20\ \mathrm{m} ;
- angles initiaux : -0{,}15^\circ, -0{,}25^\circ, -0{,}35^\circ, -0{,}50^\circ ;
- pas d’Euler : ds=0{,}05\ \mathrm{m} ;
- fenêtre horizontale : 0\leq x\leq 1100\ \mathrm{m}.
Voir le code Python
import numpy as np
import matplotlib.pyplot as plt
T_inf = 303.0
Delta_T = 20.0
H = 0.20
n_inf = 1.000265
Delta_n = (n_inf - 1) * Delta_T / T_inf
n0 = n_inf - Delta_n
def n_of_z(z):
return n_inf - Delta_n * np.exp(-z / H)
def dn_dz(z):
return (Delta_n / H) * np.exp(-z / H)
def cumulative_trapezoid(y, x):
out = np.zeros_like(x)
out[1:] = np.cumsum(0.5 * (y[1:] + y[:-1]) * np.diff(x))
return out
def turning_height(C):
"""Renvoie la hauteur de retournement si elle existe."""
if C <= n0:
return None
if C >= n_inf:
return None
ratio = (n_inf - C) / Delta_n
if ratio <= 0 or ratio >= 1:
return None
return -H * np.log(ratio)
def regularized_integrand(u, z_star, C):
"""Intégrande après le changement de variable z = z_star + u^2."""
z = z_star + u**2
value = n_of_z(z)**2 / C**2 - 1
out = np.empty_like(u)
mask = u > 0
out[mask] = 2 * u[mask] / np.sqrt(value[mask])
out[~mask] = 2 / np.sqrt(2 * dn_dz(z_star) / C)
return out
def ray_from_quadrature(z0=1.2, phi0_deg=-0.25, z_top=1.5, n_points=900):
phi0 = np.deg2rad(phi0_deg)
C = n_of_z(z0) * np.cos(phi0)
z_star = turning_height(C)
if z_star is None or z_star >= z0:
# Pas de retournement avant la route : on descend jusqu'à z=0.
z_desc = np.linspace(z0, 0.0, n_points)
value = np.maximum(n_of_z(z_desc)**2 / C**2 - 1, 0)
integrand = 1 / np.sqrt(value)
h_desc = z0 - z_desc
x_desc = cumulative_trapezoid(integrand, h_desc)
return x_desc, z_desc
# Au retournement, l'intégrande en z diverge comme 1/sqrt(z-z_star).
# Le changement de variable z = z_star + u^2 supprime cette singularité.
u_desc = np.linspace(0.0, np.sqrt(z0 - z_star), n_points)
integrand_desc = regularized_integrand(u_desc, z_star, C)
x_from_turn = cumulative_trapezoid(integrand_desc, u_desc)
x_turn = x_from_turn[-1]
z_desc = z_star + u_desc[::-1]**2
x_desc = x_turn - x_from_turn[::-1]
u_asc = np.linspace(0.0, np.sqrt(z_top - z_star), n_points)
integrand_asc = regularized_integrand(u_asc, z_star, C)
z_asc = z_star + u_asc**2
x_asc = x_turn + cumulative_trapezoid(integrand_asc, u_asc)
x = np.concatenate([x_desc, x_asc])
z = np.concatenate([z_desc, z_asc])
return x, z
def ray_from_euler(z0=1.2, phi0_deg=-0.25, ds=0.05, x_max=1100, z_top=1.5):
x = 0.0
z = z0
phi = np.deg2rad(phi0_deg)
xs = [x]
zs = [z]
n_steps = int(1.2 * x_max / ds)
for _ in range(n_steps):
n = n_of_z(z)
dx_ds = np.cos(phi)
dz_ds = np.sin(phi)
dphi_ds = (dn_dz(z) / n) * np.cos(phi)
x = x + dx_ds * ds
z = z + dz_ds * ds
phi = phi + dphi_ds * ds
xs.append(x)
zs.append(z)
if z < 0 or z > z_top or x > x_max:
break
return np.array(xs), np.array(zs)
def I_of_n(n, C):
a = n_inf
argument = ((a + C) * (n - C)) / ((a - C) * (n + C))
argument = np.clip(argument, 0, 1 - 1e-14)
return 2 / np.sqrt(a**2 - C**2) * np.arctanh(np.sqrt(argument))
def ray_from_exact_integral(z0=1.2, phi0_deg=-0.25, z_top=1.5, n_points=900):
phi0 = np.deg2rad(phi0_deg)
C = n_of_z(z0) * np.cos(phi0)
z_star = turning_height(C)
if z_star is None or z_star >= z0:
z_desc = np.linspace(z0, 0.0, n_points)
x_desc = H * C * (I_of_n(n_of_z(z0), C) - I_of_n(n_of_z(z_desc), C))
return x_desc, z_desc
z_desc = np.linspace(z0, z_star, n_points)
x_desc = H * C * (I_of_n(n_of_z(z0), C) - I_of_n(n_of_z(z_desc), C))
x_turn = H * C * I_of_n(n_of_z(z0), C)
z_asc = np.linspace(z_star, z_top, n_points)
x_asc = x_turn + H * C * I_of_n(n_of_z(z_asc), C)
x = np.concatenate([x_desc, x_asc])
z = np.concatenate([z_desc, z_asc])
return x, z
angles = [-0.15, -0.25, -0.35, -0.50]
z0 = 1.2
z_top = 1.5
x_max = 1100
ds_euler = 0.05
methods = [
("Euler pas à pas", ray_from_euler),
("Quadrature numérique", ray_from_quadrature),
("Formule analytique du modèle approché", ray_from_exact_integral),
]
fig, axes = plt.subplots(3, 1, figsize=(8, 10.8), sharex=True, sharey=True)
for ax, (title, method) in zip(axes, methods):
for angle in angles:
if method is ray_from_euler:
x, z = method(z0=z0, phi0_deg=angle, ds=ds_euler,
x_max=x_max, z_top=z_top)
else:
x, z = method(z0=z0, phi0_deg=angle, z_top=z_top)
ax.plot(x, z, label=f"{angle}°")
ax.axhline(0, linewidth=3)
ax.set_title(title)
ax.set_xlim(0, x_max)
ax.set_ylim(0, z_top)
ax.set_ylabel("Hauteur z (m)")
ax.grid(True)
axes[-1].set_xlabel("Distance horizontale x (m)")
axes[0].legend(title="Angle initial", loc="upper right")
param_text = (
f"T∞={T_inf:.0f} K ; ΔT={Delta_T:.0f} K ; H={H:.2f} m ; "
f"n∞={n_inf:.6f} ; z0={z0:.2f} m ; ds Euler={ds_euler:.2f} m"
)
fig.suptitle("Rayons au-dessus d’une route chaude", y=0.995)
fig.text(0.5, 0.965, param_text, ha="center", va="top", fontsize=9)
fig.tight_layout(rect=[0, 0, 1, 0.94])
plt.show()
La méthode d’Euler suit directement l’évolution locale du rayon, pas à pas. La quadrature numérique et la formule analytique utilisent plus directement la constante de Snell-Descartes. Les trois aboutissent toutefois pratiquement aux mêmes trajectoires.
La méthode d’Euler est plus générale : elle reste utilisable si le profil d’indice devient plus compliqué, s’il dépend aussi de x, ou si l’on ajoute des perturbations. En revanche, elle dépend du choix du pas ds et accumule une erreur numérique. Ici, elle est conservée dans sa forme classique, sans correction.
La suite dans un prochain article ! Nous exploiterons le modèle construit.
Formulaire résumé
Température :
T(z)=T_\infty+\Delta T e^{-z/H}.Indice :
n(z)=1+(n_\infty-1)\frac{T_\infty}{T(z)}.Approximation :
n(z)\simeq n_\infty-\Delta n e^{-z/H}.avec :
\Delta n=(n_\infty-1)\frac{\Delta T}{T_\infty}.Loi de Snell-Descartes continue :
n(z)\cos\varphi=C.Constante initiale :
C=n(z_0)\cos\varphi_0.Angle local :
\cos\varphi=\frac{C}{n(z)}.Pente du rayon :
\frac{dz}{dx} = \pm\sqrt{\frac{n(z)^2}{C^2}-1}.Point de retournement :
n(z_\ast)=C.