\documentclass{article}
\usepackage{a4wide}
\usepackage{graphicx}
\usepackage{multicol}
         \renewcommand{\floatpagefraction}{0.7}
         \tolerance=8000
         \pretolerance=4000




\title{Calcul débranché de conduction thermique dans un sol homogène à partir de données CMIP\author{Frédéric Hourdin}\date{\today}}


\begin{document}

\maketitle

\section{Introduction}

Ce texte présente la mise en place d'un calcul de conduction thermique dans un sol homogène forcé par des données atmosphériques. Ces données sont issues en pratique de simulations climatiques (CMIP) réalisées ici avec le modèle de l'IPSL.
Le but de l'exercice est de servir de point d'appui pour la validation des codes thermiques d'installation de surface. En effet, dans le cas particulier d'un sol homogène, on peut vérifier qu'on est capable de retrouver les températures de surface issues des simulations climatiques.

Cet exercice a permis de mettre au jour quelques difficultés à avoir en tête, à la fois dans l'extraction des données et leur utilisation.

Ce document a aussi pour vocation de guider les collègues (notamment ceux du "projet Solaire") dans l'utilisation des données issues des simulations CMIP.
Il peut également permettre d'alimenter la réflexion si on se projette dans l'utilisation de ces données pour les services climatiques.

\section{Le modèle physique}

\subsection{Les flux atmosphériques à la surface}

% Macros des constantes specifiques
\def\cs{C_s}
\def\cond{\lambda_s}
\def\sw{S}
\def\swdn{S^{\downarrow}}
\def\swup{S^{\uparrow}}
\def\lwdn{{Th}^{\downarrow}}
\def\lwup{{Th}^{\uparrow}}
\def\ts{T_s}

\def\ta{T_a}
\def\rh{RH}
\def\wpmc{W~m$^{-2}$}
\def\modv{|V_a|}
\def\cd{C_d}
\def\cp{C_p}
\def\arid{\beta}
\def\qsat#1{q_{\mbox{sat}}(#1)}
\def\qss{\qsat{\ts}}
\def\qsa{\qsat{\ta}}
\def\epot{E_{\mbox{pot}}}
\def\qa{q_a}
\def\rd{R_d}
\def\rv{R_v}
\def\tsol{T_S}
\def\fc{F_c}
\def\zn{z^*}
\def\rcdv{\kappa}
\def\rcdvmin{{\rcdv}_{min}}
\def\rcdvmax{{\rcdv}_{max}}
\def\coefev{{\kappa}_{ev}}
\def\Fatm#1{F_{a^\downarrow(#1)}}
\def\Fainsens{F_{a,0}}
\def\gama{\gamma_a}
\def\Fsurf#1{F_{s}^\downarrow(#1)}
\def\Fsinsens{F_{s,0}}
\def\gams{\gamma_s}
\def\Rdown{R^\downarrow}
\def\msol#1{d_{#1}^-}
\def\psol#1{d_{#1}^+}
\def\asol#1{A_{#1}}
\def\bsol#1{B_{#1}}
\def\dsol#1{D_{#1}}

% Macros générales 
\def\fig#1{Fig.~\ref{fg:#1}}
\def\tb#1{Table~\ref{tb:#1}}
\def\eq#1{Eq.~\ref{eq:#1}}
\def\sec#1{Section~\ref{sec:#1}}
\def\dep#1{\left(#1\right)}
\def\depb#1{\left[#1\right]}
\def\depc#1{\left\{#1\right\}}
\def\der#1#2{\frac{\partial #1}{\partial #2}}

\def\coefprem{a}
\def\tnum#1{{\tsol}_{#1}}
\def\znn#1{{\zn}_{#1}}
\def\dzn#1{{\delta \zn}_{#1}}
\def\fcn#1{{F}_{#1}^*}
\def\tslu{\ts^{lu}}
\def\tprev#1{\tnum{k}^{t-\delta t}}




On considère donc un sol homogène de capacité calorifique $\cs$ (en J~m$^{-3}$~K$^{-1}$) et de conductivité $\cond$ (en W~m$^{-1}$~K$^{-1}$). Ce sol échange avec l'atmosphère au travers de flux radiatif solaire\footnote{On sépare souvent flux "ondes courtes" et "ondes longues" plutôt que solaires et thermiques (SW et LW en anglais), mais ce qui compte en général pour la modélisation c'est de distinguer un rayonnement dont la source est extérieure (le soleil), uniquement absorbé et diffusé par l'atmosphère et la surface, et un flux thermique à la fois émis, absorbé et diffusé par l'atmosphère et la surface, les deux pouvant très bien avoir une zone de recouvrement spectral non nul.} d'ondes courtes et un rayonnement thermique ondes longues et de flux turbulents, dit sensible (noté traditionnellement $H$) et latent $LE$ (noté ainsi car produit de la chaleur latente d'évaporation par le flux d'eau $E$ associé à l'évaporation).


On notera $\swdn$, $\swup=A\swdn$ (où $A$ est l'albedo de la surface) les flux solaires descendants et montants à la surface, $\lwdn$ et $\lwup=\sigma\ts^4$ (où on suppose l'émissivité du sol égale à 1, où $\sigma=5.67\times10^{-8}$\wpmc$K^{-4}$ est la constante de Stefan Boltzmann et $\ts$ la température de surface.

Le flux sensible s'exprime sous la forme $H=\rcdv (\ts-\ta)$ où $\ta$ est la température de l'air à un niveau donné $z$ au dessus de la surface. Le coefficient d'échange s'exprime comme $\rcdv=\rho\cd\cp\modv$ où $\rho$ est la masse volumique de l'air, $\cp\simeq$1005~J~K${-1}$~kg$^{-1}$ la chaleur spécifique à pression constante de l'air, $\modv$ le module du vent horizontal au niveau $z$ et où $\cd$ est dit coefficient de trainé.
Cette formule est valable quelque soit le niveau où on considère les variables météorologiques dans les premières dizaines ou centaines de mètres typiquement, dans ce qu'on appelle couche limite de surface, définie par le fait que les flux turbulents y sont quasiment constants. Le coefficient de traînée n'a pas du coup la même valeur suivant le niveau où on considère les variables météorologiques.

Le flux évaporatif peut s'exprimer pour sa part  comme $E=\arid\epot$ où $\epot$ est l'évaporation dite potentielle, celle d'une surface d'eau libre dans les mêmes conditions d'humidité, de température et de pression, et $0 \le \arid\le 1$  est un coefficient dit d'aridité, rendant compte de la disponibilité de l'eau à la surface pour l'évaporation. L'évaporation potentielle elle même s'exprime de façon similaire au flux sensible comme
$\epot=\rho\cd\modv(\qss-\qa)$ où $\qsat{T}$ est l'humidité spécifique à saturation à la température $T$ et $\qa$ est l'humidité spécifique de l'air.
Ce qui contrôle l'évaporation, sur une surface d'eau libre, c'est le contraste d'humidité entre la couche saturée millimétrique en équilibre avec la surface d'eau libre et l'air.
Cette évaporation potentielle peut également s'approcher en introduisant l'humidité à saturation à la température de l'air et en remarquant que
   \begin{eqnarray}
   \qss-\qa&=& \qss-\qsa+\qsa-\qa\\
           &\simeq& \frac{\partial \qsat{T}}{\partial T}(\ts-\ta)-\qsa(1-\rh)\\
   	   &=& \qsa \depb{ \frac{L}{\rv \ta^2} (\ts-\ta) + 1 - \rh }
   \end{eqnarray}
On distingue par cette écriture un effet lié à la différence de température (et donc d'humidité à saturation) entre la surface et l'air au-dessus et  un effet lié à l'humidité relative de l'air.

%------------------------------------------------------------------
\subsection{Conduction thermiue  dans le sol}
%------------------------------------------------------------------

L'équation de l'évolution de la température sous la surface s'écrit
   \begin{equation}
   \cs\der{\tsol}{t}=-\der{\fc}{z}
   \end{equation}
où $z$ est la profondeur mesurée depuis la surface vers le fond et où le flux conductif s'écrit
   \begin{equation}
   \fc=-\cond\der{\tsol}{z}
   \end{equation}
Sous la surface, l'équation d'évolution de la température peut se simplifier en introduisant le changement de variable $\zn=z\sqrt{\cs/\cond}$:
   \begin{equation}\label{eq:normalisee}
   \der{\tsol}{t}=\frac{\cond}{\cs}\der{^2\tsol}{z^2}=\der{^2\tsol}{{\zn}^2}
   \end{equation}
   Avec cette coordonnée homogène à la racine carrée d'un temps, les flux d'énergie s'expriment comme 
   \begin{equation}
   \fc=-I\der{\tsol}{\zn}
   \end{equation}
   où $I=\sqrt{\cond\cs}$ appelée inertie thermique est le seul paramètre entrant en compte dans le fonctionnement de la température de surface, vue depuis l'atmosphère. Le comportement de la propagation et de l'atténuation de ondes thermiques sous la surface est indépendant des propriétés du sol dans la coordonnée $\zn$ mais la connaissance indépendante de $\cond$ et $\cs$ est indispensable pour attribuer une vraie profondeur physique à la température estimée à un niveau $\zn$.

Les valeurs typiques des conductivités et des capacités calorifique du sol tournent autour de $\cond\simeq$1~W~m$^{-1}$~K$^{-1}$ et $\cs\simeq 1\times 10^6$~J~m$^{-3}$~K$^{-1}$, correspondant donc à une inertie thermique $I\simeq$1000~J~m$^{-2}$~s$^{-1/2}$~K.
Les valeurs typiques pour les différents types de sol peuvent être trouvées sur~:
\begin{verbatim}
https://energieplus-lesite.be/donnees/enveloppe44/caracteristiques-thermiques-des-sols/
\end{verbatim}
	

%-------------------------------------------------------------------------------
\subsection{Continuité des flux à la surface\label{sc:surface}}
%-------------------------------------------------------------------------------

La continuité des flux en surface, orientés vers le bas comme le flux conductif sous la surface, s'écrit finalement
\begin{eqnarray}
-I\der{\tsol}{\zn}&=&\swdn-\swup+\lwdn-\lwup-H-LE\label{eq:flux1}\\
		  &=&(1-A)\swdn+\lwdn-\sigma{\ts}^4-\rcdv(\ts-\ta)-\beta\rcdv\frac{L}{\cp}(\qss-\qa)\label{eq:flux2}\\
		  &\simeq& \Rdown -\sigma{\ts}^4-(\rcdv+\coefev)\times(\ts-\ta)-LE_0\label{eq:decompE}
\end{eqnarray}
avec
\begin{eqnarray}
\Rdown&=&(1-A)\swdn+\lwdn\\
LE_0&=&\beta\rcdv\frac{L}{\cp}\qsa(1-\rh)\\
\coefev&=&\beta\rcdv\frac{L^2}{R {\ta}^2\cp}\qsa
\end{eqnarray}
Le passage de l'\eq{flux1} à l'\eq{flux2} part du principe  que le coefficient de trainée est le même pour la température et l'humidité (mais on peut aussi imaginer un sens élargi de $\beta$ si ce n'est pas le cas) de sorte que
$\L\rho\cd\modv=\rcdv\L/\cp$. L'\eq{decompE} permet de faire apparaître toutes les dépendances en température de surface si on suppose que la température et l'humidité de l'air à la surface sont connues, données par les simulations climatiques. Cette approximation tombera si l'installation thermique au sol a une emprunte spatiale suffisante pour que le passage de l'air sur l'installation lui laisse le temps de s'ajuster localement. Ce cas est au c{\oe}ur du travail de thèse de Louis d'Alençon et de du projet MC2 d'étude du couplage ville climat financé par l'Agence Nationale de la Recherche.

Mais ici, où on essaie de vérifier les différents modèles utilisés pour retrouver la température de surface simulée par le modèle, on supposera ces variables météorologiques connues.



%===============================================================================
\section{Les données météorologiques}
%===============================================================================

%-------------------------------------------------------------------------------
\subsection{Les simulations CMIP}
%-------------------------------------------------------------------------------

On va utiliser ici des simulations réalisées dans le cadre du Projet d'Intercomparaison des Modèles de Couplés de climat (CMIP) en amont des rapports du Giec (IPCC en anglais).
Ces simulations sont orchestrées au niveau international avec des protocoles communs très précis et les données de sorties des simulations sont mises à disposition avec une nomenclature très précise de d'arborescence de fichier commune à tous les centres.

Les données des simulations réalisées par la plupart des groupes sont libres et distribuées au travers de l'ESGF, une infrastructure développée pour l'occasion.

On va considérer plus spécifiquement les simulations réalisées avec le modèle de climat de l'IPSL, dans la mesure où on va chercher à utiliser des sorties à très haute fréquence temporelle permettant de reconstituer l'évolution des températures depuis le cycle diurne jusqu'aux évolutions centenales (notamment en réponse à l'augmentation des gaz à effet de serre).


%-------------------------------------------------------------------------------
\subsection{Accès aux données}
%-------------------------------------------------------------------------------


Les variables sont disponibles sur le cluster ciclad de l'IPSL, sur le répertoire
\begin{verbatim}
/bdd/CMIP6/[sous-projet]/[institut]/[configuration]/[experience]/[realisation]/[frequence]
\end{verbatim}

Pour les simulations utilisées ici, on prend des simulations réalisées avec le modèle de l'IPSL:\\
\begin{verbatim}
institiut=IPSL
confiugration=IPSL-CM6A-LR
\end{verbatim}


\begin{table}
\centerline{\begin{tabular}{lll}\hline
CMIP & piControl & Simulation de contrôle ; forçages pré-industriels  \\
     & historical & Simulations historiques, reconstructions 1850-2014\\
     & historical-EXT & Simulations historiques étendues à 2060 (IPSL)  \\
     & 1ptCO2 & Simulations avec augmentation de 1$\%$ du CO$_2$ par an\\
     & abrupt-4xCO2 & Multiplication instantanée par 4 du CO$_2$ \\
     & amip & Reconstructions avec atmosphère seule \\\hline
DAMIP & hist-aer & Reconstructions ne tenant compte que du changement des aérosols \\
      & hist-GHG & Reconstructions ne tenant compte que des gaz à effet de serre  \\
      & hist-nat & Reconstructions ne tenant compte que des "Forçages naturels"" \\\hline
ScenarioMIP & ssp126       &    \\
            & ssp370       &    \\
            & ssp119       &    \\
            & ssp126       &    \\
            & ssp245       &    \\
            & ssp370       &    \\
            & ssp434       &    \\
            & ssp434       &    \\
            & ssp460       &    \\
            & ssp460       &    \\
            & ssp534-over       &    \\
            & ssp585       &    \\
\end{tabular}}
\caption{Configurations disponibles avec des sorties tri-horaires\label{tb:configs}}
\end{table}

Les sous-projets et configurations possibles possibles sont listés dans la \tb{configs}.
Pour chaque simulation, un centre peut fournir un "ensemble" de simulations.
Chaque membre de cet ensemble diffère généralement par l'utilisation d'un état initial différent. Chaque "relase" est baptisée en générale
\begin{verbatim}
realisation=rXi1p1f*
\end{verbatim}
Ici, on utilise la plus haute fréquence temporelle disponible, à savoir
\begin{verbatim}
frequence=3hr
\end{verbatim}


%-------------------------------------------------------------------------------
\subsection{Extraction des données}
%-------------------------------------------------------------------------------

\begin{table}
\centerline{\begin{tabular}{llll}
CMIP & Texte & description & unité \\\hline
tslsi  & $\ts$ & Surface temperature & K\\
tas  & $\ta$ & Air 2m temperature & K \\
huss  & $\qa$ & Air 2m specific humidity & kg/kg \\
ps  & $P_s$  & Surface pressure & Pa \\
hfls  & $LE$ & Latente heat flux & \wpmc \\
hfss  & $H$ & Sensible heat flux & \wpmc \\
rlds  & $\lwdn$ & Downward thermal radiative flux (saied LW) & \wpmc \\
rsds  & $\swdn$ & Downward solar flux (SW) & \wpmc \\
rsus  & $\swup$ & Upward solar flux & \wpmc \\
pr  & P & Precipitation & mm/day \\
snc  & f$_{sn}$ & Couverture neigeuse & fraction (freq jour)
\end{tabular}}
\caption{Nom des variables CMIP utilisées pour ce travail. A noter que les variables météorologiques (tslsi, tas, huss, ps) sont des variables instantanées à la fin des intervalles de 3 heures, alors que les autres, de type flux, sont calculées en moyenne sur les 3 heures. La fraction neigeuse n'est elle disponible qu'à la fréquence journalière.\label{tb:nomvarmip}}
\end{table}

L'extraction des données se fait en plusieurs étapes.
D'abord sur ciclad, on extrait un point pour toutes les variables listées dans le \tb{nomvarmip} à l'aide de commandes 
\begin{verbatim}
ncks -d lon,3. -d lat,48. 
\end{verbatim}
et
\begin{verbatim}
ncrcat
\end{verbatim}
pour concaténer 
les fichiers d'origine, souvent séparés en plusieurs parties correspondant à des périodes successives.

Les variables sont ensuite mise en forme dans un fichier colonne ascii comme illustré sur la \fig{ascii}.
\begin{figure}
{\small
\begin{verbatim}
Date        Heure    Ts(K)  Ta(K) Qa(g/kg) Ps(Pa) P(mm/j) SWdn   SWup  LWdn        H      LE SnowC
01-JAN-1850 01:30:00 267.94 268.45    2.66  99752   0    0.00    0.00  235.13   -9.05   -6.15 13.01
01-JAN-1850 04:30:00 268.12 268.00    2.56  99820   0    0.00    0.00  254.72   -1.47   -3.16 13.01
01-JAN-1850 07:30:00 268.63 267.67    2.41  99939   0   25.26    5.45  256.84    8.32    2.20 13.01
01-JAN-1850 10:30:00 273.15 271.47    3.23  99990   0  230.63   48.88  265.30   39.70   23.30 13.01
01-JAN-1850 13:30:00 273.15 272.02    3.45 100060   0  163.65   34.72  288.62   30.69   19.39 13.01
01-JAN-1850 16:30:00 271.66 271.04    3.19 100195   0    3.17    0.69  305.22   11.07    5.82 13.01
01-JAN-1850 19:30:00 271.08 270.63    3.05 100280   0    0.00    0.00  296.45    7.64    4.03 13.01
01-JAN-1850 22:30:00 270.84 270.54    2.93 100361   1    0.00    0.00  295.96    6.32    4.58 13.01
02-JAN-1850 01:30:00 270.55 270.22    2.84 100390   0    0.00    0.00  293.31    4.40    5.89 18.07
02-JAN-1850 04:30:00 270.60 270.41    2.91 100427   1    0.00    0.00  292.46    3.90    5.49 18.07
02-JAN-1850 07:30:00 271.36 270.65    2.91 100506   1   19.80    4.77  293.30    4.35    6.19 18.07
02-JAN-1850 10:30:00 273.15 272.28    3.28 100481   1  191.69   46.23  297.31   28.26   23.88 18.07
02-JAN-1850 13:30:00 273.15 272.68    3.30 100456   1  131.19   31.29  304.32   15.70   17.19 18.07
02-JAN-1850 16:30:00 272.43 272.18    3.24 100453   1    3.24    0.79  304.20    2.85    7.56 18.07
02-JAN-1850 19:30:00 271.77 271.80    3.29 100418   0    0.00    0.00  290.84    0.14    2.40 18.07
02-JAN-1850 22:30:00 271.09 271.27    3.19 100327   0    0.00    0.00  282.45   -0.92    1.15 18.07
03-JAN-1850 01:30:00 271.22 271.16    3.16 100203   0    0.00    0.00  284.72   -0.64    0.72 30.35
03-JAN-1850 04:30:00 271.34 271.25    3.13 100080   0    0.00    0.00  295.25    1.09    3.21 30.35
\end{verbatim}
}
\caption{Organisation en colonne du fichier météo correspondant à une simulation de l'archive CMIP\label{fg:ascii}}
\end{figure}

%-------------------------------------------------------------------------------
\subsection{Calcul de $\rcdv$ et $\arid$}
%-------------------------------------------------------------------------------

La variable $\rcdv$ n'est pas disponible dans les archives CMIP.
Même si elle l'était, celle du modèle est calculée par rapport à l'altitude  de la première couche. Ici, on choisit d'utiliser les variables "à 2m" (niveaux standard pour les observation et les sorties des températures et humidité de l'air en surface, alors que le vent est "sorti" à 10m) et de recalculer le coefficient à partir des autres variables $\rcdv=H/(\ts-\ta)$.

C'est simple a priori.
Dans la pratique on voit bien que c'est problématique déjà à cause de la division par la différence de température qui peut être nulle.
En plus, les températures correspondent à des valeurs instantanée
alors que le flux $H$ est moyenné sur 3 heures entre deux instants.

En pratique, on impose que le $\rcdvmin <\rcdv<\rcdvmax$ en \wpmc/K.  Les valeurs $\rcdvmin=0.2$ et  $\rcdvmax=50$ ont été choisies après un certain nombre de tests.

En pratique, pour calculer $\rcdv$, on commence par regarder si $H$ et $\ts-\ta$ sont du même signe quand on considère les températures moyennes en début et en fin  de pas de temps
Si ce n'est pas le cas, on regarde c'est le cas avec les températures prises soit en début ou en fin de pas de temps.
Une fois ainsi calculer la valeur de $\rcdv$, on la borne.

De la même façon, le calcul de $\arid$ est calculé comme 
$\arid=E/\epot$ avec $\epot=\rcdv L/\cp(\qss-\qa)$.


\begin{figure}
\centerline{\includegraphics[angle=-90,width=13cm]{France1.1700.kz0.2.50.f8.d113100920.pdf}}
\centerline{\includegraphics[angle=-90,width=13cm]{France1.1700.kz0.2.50.f8.d113100921.pdf}}
\caption{Validation du calcul des coefficients d'échange en surface sur deux années consécutives. En haut, comparaison du flux sensible $H$ issu de la simulation (courbe noire épaisse) et du flux recalculé comme le produit $\rcdv(\ts-\ta)$ (rouge) où le coefficient $\rcdv$ est lui même calculé comme $H/(\ts-\ta)$ selon la procédure expliquée dans le texte. En bas, le flux $LE$ issu de la simulation CMIP (noir) est comparé au flux recalculé comme $\beta\rcdv (L/\cp)(\qss-\qa)$ (rouge). Le coefficient $100\times\beta$ est montré en cyan. Le flux en vert est la partie liée à $\rh$ de la décomposition de l'\eq{decompE}
$E0=\beta\rcdv\frac{L}{\cp}\qsa(1-\rh)$.
La courbe bleue marine montre le même flux auquel on a ajouté le terme $\beta\rcdv\frac{L}{\cp}\qsa\frac{L}{\rv (\ta)^2}(\ts-\ta)$.
Pour aider à interpréter les variations de $\beta$ on montre également la précipitation en mm/jour (magenta). Mis à part les précipitations montrées à une fréquence horaire, toutes les variables sont montrées en moyenne glissante sur 5 jours. \label{fg:meteo1}}
\end{figure}


\begin{figure}
\centerline{\includegraphics[angle=-90,width=13cm]{France1.1700.kz0.2.50.f8.d113100930.pdf}}
\centerline{\includegraphics[angle=-90,width=13cm]{France1.1700.kz0.2.50.f8.d113100931.pdf}}
\caption{Même figure que la \fig{meteo1} mais sur une période plus courte et sans lissage en temps, pour vérifier la reconstruction de la haute fréquence, et illustrer l'importance du cycle diurne sur ces flux. \label{fg:meteo2}}
\end{figure}

Avant d'entrer dans le vif des intégrations numériques, on montre sur les \fig{meteo1} avec un lissage à 5 jours, et la \fig{meteo2} sans lissage (fréquence 3 heures) la comparaison des flux $H$ et $LE$ et issus des simulations CMIP et reconstitués en fonction des variables $\ta$, $\ts$ et $\qa$ après inversion des coefficient $\rcdv$ et $\beta$.

On voit quelques oscillations pathologiques du coefficient $\rcdv$ qui n'empêchent pas de reconstituer correctement les flux.
On voit aussi les effets de compensation entre flux sensible et latent entre les deux années sur la \fig{meteo1} avec une année 1857 plus sèche que l'année 1858 l'été pour cette simulation particulière.

%===============================================================================
\section{Intégration numérique en différences finies} 
%===============================================================================

\subsection{Discrétisation spatiale}

\begin{figure}
\begin{verbatim}
    !
    !   ~ / ~ / ~ / ~ / ~ / ~ / ~ / ~ / ~   Atmosphère 
    !   ---------------------------------   Surface Ts
    !
    !   + + + + + + + + + + + + + + + + +    Premeière couche  T(1), dz(1)
    !                                           
    !   ---------------------------------   
    !
    !   + + + + + + + + + + + + + + + + + \  Milieu de couche T(k-1), dz(k-1)
    !                                      \
    !   ---------------------------------    Flux Fc(k) entre les couches k-1 et k
    !                                        (Epaisseur [dz(k-1)+dz(k)]/2
    !                                      /
    !   + + + + + + + + + + + + + + + + + /  Milieu de couche T(k), dz(k)
    !
    !
    !   ---------------------------------
\end{verbatim}
\caption{Schéma de la répartition des niveaux verticaux et des variables dans le modèle discrétisé de sol.\label{fg:discret}}
\end{figure}

On discrétise le sol en $N$ couches d'épaisseur variable $\dzn{k}$.
En pratique on prend $\dzn{k+1}=\xi\dzn{k}$, avec $\xi=1.3$.
Les profondeur des interfaces entre couches sont donc données par 
$\znn{1}=0.$, $\znn{k+1}=\Sigma_{k'=1}^n \dzn{k'}$.
A ces interfaces, on va calculer les flux conductifs.
Comme indiqué sur la \fig{discret}, les températures seront elles calculées au milieu des couches (ou en moyenne sur la couche), de sorte que la distance entre deux milieux de couches (entre deux niveaux de température) vaut $[\dzn{k-1}+\dzn{k}]/2$.


Avec ces notations, les flux conductifs s'écrivent :
\begin{equation}\label{eq:fluxnum}
\frac{\fc}{I}(\znn{k-1/2})=\fcn{k}=-\frac{\tnum{k}-\tnum{k-1}}{[\dzn{k}+\dzn{k-1}]/2}
\end{equation}

Le laplacien de l'\eq{normalisee} peut alors se réécrire
\begin{equation}\label{eq:evolnum}
\der{^2\tsol}{{\zn}^2}(k)\simeq
- \frac{ \fcn{k+1} - \fcn{k} }{ \dzn{k} }
\end{equation}

\begin{table}
\centerline{\begin{tabular}{r|r|r|r}
       & \multicolumn{2}{|c|}{Profondeur $\zn$ (s$^{-1/2}$)} & Période de profondeur  \\
       & \multicolumn{2}{|c|}{$\simeq z$ en mm} & de pénétration \\
Niveau & interfaces & milieux         & $z$, en jour \\\hline
  1 &      0.00 &     15.00 &       0.02 \\
  2 &     30.00 &     49.50 &       0.18 \\
  3 &     69.00 &     94.35 &       0.65 \\
  4 &    119.70 &    152.65 &       1.69 \\
  5 &    185.61 &    228.45 &       3.79 \\
  6 &    271.29 &    326.99 &       7.77 \\
  7 &    382.68 &    455.08 &      15.06 \\
  8 &    527.49 &    621.61 &      28.10 \\
  9 &    715.73 &    838.09 &      51.07 \\
 10 &    960.45 &   1119.52 &      91.13 \\
 11 &   1278.58 &   1485.37 &     160.43 \\
 12 &   1692.16 &   1960.98 &     279.61 \\
 13 &   2229.81 &   2579.28 &     483.73 \\
 14 &   2928.75 &   3383.06 &     832.20 \\
 15 &   3837.38 &   4427.98 &    1425.67 \\
 16 &   5018.59 &   5786.38 &    2434.55 \\
 17 &   6554.16 &   7552.29 &    4147.28 \\
 18 &   8550.41 &   9847.97 &    7051.80 \\
 19 &  11145.54 &  12832.37 &   11973.46 \\
 20 &  14519.20 &  16712.07 &   20308.00
\end{tabular}}
\caption{Exemple de discrétisation d'un modèle de sol
pour $\cond\simeq$1~W~m$^{-1}$~K$^{-1}$ et $\cs\simeq 1\times 10^6$~J~m$^{-3}$~K$^{-1}$.\label{tb:zsol}}
\end{table}

Pour les calculs en différence finie décrit plus bas, 
Avec une première couche à 30~s$^{-1}$, $\xi=1.3$ et $N=20$, la première couche du modèle correspond à un temps de 700~s, les 4 premières couches sont dans la profondeur de pénétration du cycle diurne, la couche 14 correspond à la pénétration du cycle annuel et la couche 20 à une 30 aines d'années, soit une profondeur de sol de typiquement 10~m.

En calculant les solutions de l'équation de conduction pour une condition en flux oscillant en surface (cf. \fig{analytique}), on peut montrer que la profondeur de pénétration $z_P$ d'une onde de période $P$ (profondeur à laquelle sont amplitude est réduite d'un facteur $e$) est donnée par
\begin{equation}
z_P =\sqrt{\frac{\cond P}{2\cs \pi}}
\end{equation}

\begin{figure}
\centerline{\includegraphics[width=14cm]{DocScan_08_29_2021-crop.pdf}}
\caption{Calcul des solution oscillantes de la température de surface pour un flux en surface sinusoïdal imposé.\label{fg:analytique}}
\end{figure}

Exercice : le facteur 2 manque, page 39 de ma thèse. Vérifier qui a raison.

On voit donc que la coordonnée normalisée correspondant à la pénétration d'une onde de période $P$ est ${\zn}_P=\sqrt{P/(2\pi)}$.
Dans ces coordonnées, la profondeur de pénétration de l'onde diurne est typiquement de $\simeq 115$~s$^{-1}$, soit une profondeur de 0.115~m pour $\cond\simeq$1~W~m$^{-1}$~K$^{-1}$ et $\cs\simeq 1\times 10^6$~J~m$^{-3}$~K$^{-1}$.



\subsection{Intégration temporelle}

On peut écrire l'intégration temporelle en temps en introduisant une notation pour l'évolution des variables sur un pas de temps $\Delta_t$
\begin{equation}
\Delta_t\tnum{k} = ({\fcn{k}} - {\fcn{k+1}}) \frac{\delta t }{ \dzn{k} }
\end{equation}
Dans un schéma 
explicite en temps, plus "intuitif" ou "direct", mais facilement instable pour ce type d'équations, on estime le membre de droite au début du pas de temps.
Dans un schéma implicite, inconditionnellement stable pour cette équation, mais nécessitnat l'inversion d'une matrice (heureusement tridiagonale), on estime le membre de droite à la fin du pas de temps.
Si on suppose que $t$ est l'intstant où on évalue le terme de droite de l'équation, la distinction revient à choisir entre les deux définitions de $\Delta_t$:
\begin{eqnarray}
\Delta_t \tnum{k} &=& \tnum{k}^{t+\delta t} - \tnum{k}^{t}\mbox{    schéma explicite}\\
\Delta_t \tnum{k}  &=& \tnum{k}^{t} - \tnum{k}^{t-\delta t}\mbox{    schéma implicite}
\end{eqnarray}
le flux lui même s'écrivant comme:
\begin{equation}
\fcn{k}=\frac{\tnum{k-1}^t-\tnum{k}^t}{[\dzn{k}+\dzn{k-1}]/2}
\end{equation}
pour les interfaces entre couches indicées de $2$ à $N$.
On supposera que le flux au fond du modèle de sol est nul
\begin{equation}
\fcn{N+1}=0
\end{equation}

\subsection{Continuité des flux à la surface}

A la surface, où on écrit la continuité avec les flux atmosphériques, on peut faire les trois hypothèses suivantes~:

\underline{Hyp1:} on considère que la température de surface vue par l'atmosphère est directement la température de la première couche.
A ce moment là, on impose directement le flux atmosphérique en haut de la première couche
\begin{equation}
I \fcn{1}=\Fatm{\tnum{1}}
\end{equation}

\underline{Hyp2:} on peut être tenté (je le dis comme ça parce que c'est la première version que j'ai codée, et que ça ne parait pas si évident) d'interpoler
la température à la surface linéairement à partir des température des première et seconde couche, ce qui s'écrit ici
\begin{equation}\label{eq:fluxnum}
\ts^{ext}=\tnum{1}+\alpha\dep{\tnum{1}-\tnum{2}} \mbox{  avec  }\alpha=\frac{\dzn{1}}{\dzn{1}+\dzn{2}}
\end{equation}
\begin{equation}
I \fcn{1}=\Fatm{\ts^{ext}}
\end{equation}
et d'écrire à nouveau la continuité des flux en surface.
Mais cette interpolation linéaire, dans un sol homogène, dit fondamentalement  que le flux conductif
est le même aux interfaces haute et basse de la première couche du modèle, et donc fondamentalement que la température de cette couche ne peut pas évoluer.
C'est sans doute à la fois vrai pour une couche suffisamment fine au sens où (et on dit la même chose fondamentalement dans la couche limite atmosphérique) une différence "significative de flux (par exemple si on mettait un flux nul en bas de la première couche) aboutirait à une évolution temporelle délirante de la température de cette couche.
En fait, l'évolution de la température d'une telle couche est beaucoup plus lente que ce que ne pourrait laisser présager l'ordre de grandeur des flux parce que ces flux se répartissent presque immédiatement dans une couche plus épaisse.

\underline{Hyp3:} on considère que la température de surface est différente de la température au centre de la première couche et on peut alors écrire
le flux en haut de la première couche comme
\begin{equation}
\fcn{1}=\frac{\ts^t-\tnum{1}^t}{\dzn{1}/2}
\end{equation}
où $\ts$ est donné par
\begin{equation}
I \fcn{1}=\Fatm{\ts}
\end{equation}

\subsection{Linéarisation du flux atmosphérique}

Dans le cas où on doit inverser une équation pour trouver la température de surface, il est possible de linéariser le flux en surface, le plus simple consistant à linéariser par rapport au pas de temps précédent.
Etant donné qu'on connaît tous les flux dans la simulation climatique, on peut aussi jouer à prendre en compte tout ou partie des dépendances en température des différents flux.

On peut à  l'extrême imposer tous les flux, y compris l'émission thermique de la surface
\begin{equation}
\Fatm{\ts}=-\sigma {\tslu}^4 + R^\downarrow -H -LE
\end{equation}
en lisant directement $\tslu$, $R^\downarrow$, $H$ et $LE$ dans les archives météo (ici les sorties des simulations CMIP).
On peut à l'autre extrême utiliser les expressions dérivées plus haut pour $H$ et $LE$ et linéariser la fonction de Planck de sorte que 
\begin{equation}
\Fatm{\ts}=-\sigma {\ts^t}^4-4\sigma {\ts^{t-\delta t}}^3 \dep{\ts^{t}-\ts^{t-\delta t}} + R^\downarrow + \dep{\rcdv+\coefev}\times \dep{\ta - \ts^{t}} - LE_0
\end{equation}

De façon générale, on écrira
\begin{equation}
\Fatm{\ts}=\Fainsens+\gama\ts
\end{equation}
et on écrira l'égalité avec le flux sous la surface
\begin{equation}
\Fsurf{\ts}=\Fsinsens+\gams\ts
\end{equation}
d'où
\begin{equation}
\ts=\frac{\Fsinsens-\Fainsens}{\gama-\gams}
\end{equation}

\subsection{Tests du schéma explicite}

\begin{figure}
% series.sh suivi de inertie_post.sh
\centerline{\includegraphics[angle=-90,width=13cm]{fig_Inertie0.pdf}}
\centerline{\includegraphics[angle=-90,width=13cm]{fig_Inertie1.pdf}}
\caption{Effet de l'inertie thermique sur les température pour un point Saharien (0W,24N). En haut, températures en °C pour 6 jours.
En bas, amplitude (quadratique moyenne) du cycle diurne après application d'une moyenne glissante avec une fenêtre de 30 jours.\label{fg:inertie}}
\end{figure}


\begin{figure}
% series.sh suivi de inertie_post.sh
\centerline{\includegraphics[angle=-90,width=13cm]{Sahara1.1450.kz0.2.50.f2.d301900011.pdf}}
\centerline{\includegraphics[angle=-90,width=13cm]{fig_explicit1.pdf}}
\centerline{\includegraphics[angle=-90,width=13cm]{fig_explicit2.pdf}}
\caption{Tests du schémas explicite au Sahara : le panneau  du haut montre les erreurs sur la température de surface avec le schéma explicite et la grille la plus fine. Les deux panneaux  suivants montrent l'amplitude du cycle diurne lissée à un mois et la différence entre température de la surface et température de l'atmosphère, également lissée à un mois pour 5 choix numériques avec les 3 grilles et les Hyp1 et Hyp2 pour la condition de couplage à la surface.
\label{fg:explicit}}
\end{figure}

\begin{table}
\centerline{\begin{tabular}{ccccc}
$N$ & Hyp & $\delta t$ & Erreur RMS (K) & Erreur Max (K) \\\hline
\multicolumn{5}{c}{Schéma temporel explicite} \\\hline
300 & 1 &  1.2s &  1.431 &  5.782  \\
 70 & 1 &   10s &  1.448 &  5.783  \\ 
 20 & 2 &   10s &  1.418 &  5.784  \\ 
  - & 1 & 6min  &  1.576 &  6.152  \\ 
  - & 2 & 6min  &  1.391 &  5.830  \\\hline
\multicolumn{5}{c}{Schéma temporel implicite} \\\hline
 70 & 3 &  10s  &  1.420 &  5.784  \\ 
  - & - & 10min &  1.469 &  5.814  \\ 
  - & - & 20min &  1.521 &  5.943  \\ 
  - & - &  1hr  &  1.750 &  6.467  \\ 
  - & - &  3hr  &  2.503 &  8.078  \\ 
 20 & - & 6min  &  1.433 &  5.807  \\ 
  - & - & 10min &  1.453 &  5.809  \\ 
  - & - & 20min &  1.504 &  5.894  \\ 
  - & - &   1hr &  1.727 &  6.413  \\ 
  - & - &   3hr &  2.454 &  7.975  \\\hline
\end{tabular}}
\caption{Erreurs commises avec les différents schémas d'intégrations temporelles, grilles verticales et pas de temps, en prenant les températures de surface de l'archive CMIP comme référence.\label{tb:ExpImpl}}
\end{table}


On commence par présenter des résultats avec une résolution temporelle explicite.
En plus de la grille à 20 niveaux présentée plus haut, on effectue des tests avec deux grilles plus fines. 
Comme on l'a dit plus haut, le schéma explicite est instable pour des pas de temps trop grands. En pratique, le pas de temps maximum est en $(\delta z_{min})^2$.
Les trois grilles ne peuvent donc tourner qu'avec des pas de temps relativement courts.
On donne ci dessous les paramètres contrôlant les trois grilles ainsi que le pas de temps $\delta t$ retenu:
\begin{description}
\item{grille 1:}  $\dzn{1}$=30~s$^{-1}$, $\xi=1.3$, $N=20$, $\delta t$=6~min (instable jusque env. 15~min)\\ 
% niter = 30 (niter=12 instable)
\item{grille 2:}  $\dzn{1}$=5~s$^{-1}$, $\xi=1.1$, $N=70$, $\delta t$=10~s (instable jusqu'à env. 22~s)\\
% niter = 30 (niter=(500 instable, 600 OK)
\item{grille 3:}  $\dzn{1}$=2~s$^{-1}$, $\xi=1.02$, $N=300$, $\delta t$=1.2~s (instable jusqu'à env. 2.2~s)
% niter = 30 (niter=12 instable)
\end{description}

Exercice : vérifier la qu'on retrouve bien ces pas de temps critiques par l'analyse numérique.

\subsection{Inversion du schéma implicite}

\begin{figure}
% series.sh suivi de inertie_post.sh
\centerline{\includegraphics[angle=-90,width=13cm]{fig_ImplDis11.pdf}}
\centerline{\includegraphics[angle=-90,width=13cm]{fig_ImplDis12.pdf}}
\caption{Tests du schémas implicite avec la grille 1 pour différents pas de temps.
\label{fg:ImplDis1}}
\end{figure}


En combinant les équations données plus haut, on peut réécrire le schéma implicie comme 
\begin{equation}\label{eq:impl}
\tnum{k}-\tprev{k}=
\psol{k}\tnum{k+1}-\dep{\psol{k}+\msol{k}}\tnum{k}+\msol{k}\tnum{k-1}
\end{equation}
avec
\begin{eqnarray}
\msol{k}&=&\frac{\delta t}{\dzn{k}}\times\frac{2}{\dzn{k-1}+\dzn{k}}\\
\psol{k}&=&\frac{\delta t}{\dzn{k}}\times\frac{2}{\dzn{k}+\dzn{k+1}}
\end{eqnarray}

Pour la couche du fond, on a la même équation avec $\psol{N}=0$.

Pour la couche de surface, on va considérer ici l'Hyp3: on calcule le flux au bord par différence finie entre le milieu de la couche et la surface. On peut utiliser l'\eq{impl} avec
$\msol{1}=2\delta t/\dzn{1}^2$ et $\tnum{0}=\ts$. 

On va inverser la matrice par une méthode de type pivot de Gauss.
On va commencer à chercher des coefficient $\asol{k}$ et $\bsol{k}$ tels que
\begin{equation}\label{eq:recur1}
\tnum{k}=\asol{k}+\bsol{k}\tnum{k-1}
\end{equation}
(avec $\tnum{1}=\asol{1}+\bsol{1}\ts$).
On a dans la dernière couche directement
\begin{eqnarray}
\asol{N}&=& \tprev{N}/\dsol{N}\\
\bsol{N}&=& \msol{N}/\dsol{N}\\
\dsol{N}&=& 1+\msol{N}
\end{eqnarray}
et dans les suivantes, après substitution de $\tnum{k+1}$ par $\asol{k+1}+\bsol{k+1}\tnum{k}$ dans l'\eq{impl}:
\begin{eqnarray}
\asol{k}&=& \dep{\tprev{k}+\psol{k}\asol{k}}/\dsol{k}\\
\bsol{k}&=& \msol{k}/\dsol{k}\\
\dsol{k}&=& 1+\psol{k}\dep{1-\bsol{k}}+\msol{k}
\end{eqnarray}

On commence par calculer les coefficients $\asol{k}$ et $\bsol{k}$ depuis la couche $N$ jusqu'à la couche 1, puis on calcule la température de surface à partir de la continuité des flux, puis on calcule les nouvelles températures jusqu'au fond grâce à l'\eq{recur1}.
La \fig{ImplDis1} montre la simulation de la différence $\ts-\ta$ et l'amplitude du cycle diurne de $\ts$ avec ce schéma.
Les erreurs quadratiques moyennes et max sont reportés dans lee \tb{ExpImpl}.

\subsection{Prise en compte de l'évaporation}

\begin{figure}
% series.sh suivi de inertie_post.sh
\centerline{\includegraphics[angle=-90,width=13cm]{Evap.pdf}}
\centerline{\includegraphics[angle=-90,width=13cm]{fig_Evap2.pdf}}
\caption{Test des différentes façons de prendre en compte l'évaporation pour un point situé au Sahel
\label{fg:Evap}}
\end{figure}

On avait volontairement choisi un point au Sahara pour éviter de considérer la prise en compte de l'évaporation.
Après un peu de tâtonnement, on en est arrivé à la formulation proposée plus haut. Pour regarder l'effet de la prise en compte de cette évaporation, on se concentre sur 
un point au Sahel. Le Sahel est intéressant car la pluie est concentrée sur une saison courte l'été (saison des pluies ou saison de mousson).
Du coup, l'évaporation n'est significative qu'à cette époque.
La \fig{Evap} montre la capacité à la fois à reproduire les flux et leur dépendance à la température une fois inversés les coefficients $\rcdv$ et $\coefev$ et la capacité à reproduire les températures en tenant compte de ces flux évaporatif (panneau du bas). En terme d'erreur RMS (non montrée), le meilleur résultat (mais de peu) est celui où on impose l'évaporation lue (courbe rouge). Mais celle où on décompose le flux latent avec sa dépendance à la température marche bien aussi (courbe bleue). On illustre également sur cette figure ce que serait les températures si on ne tenait compte que des flux radiatifs (courbe magenta), que du flux sensible mais pas de l'évaporation (courbe verte) et enfin de la partie non dépendante à la température de l'évaporation.

Le contrôle de la température de surface par l'évaporation est évidemment important. Il est en fait largement sous estimé dans les calculs effectués ici du fait que la température atmosphérique est affectée. Pendre en compte la dépendance à l'évaporation de la température de l'air aurait un impact très fort, pas forcément sur $\ts-\ta$ mais sur la température elle-même.

On voit aussi que la décomposition du flux latent en fonction de la température fonctionne relativement bien, et qu'environ la moitié du flux peut être attribué au contraste de température entre la surface et l'air, dans cette décomposition entre effet de température et effet de l'humidité relative.



\subsection{Prise en compte de la neige}

\begin{figure}
% series.sh suivi de inertie_post.sh
\centerline{\includegraphics[angle=-90,width=13cm]{fig_Neige0.pdf}}
\centerline{\includegraphics[angle=-90,width=13cm]{fig_Neige2.pdf}}
\caption{Test d'une astuce pour traiter les jours avec neige, sur un point dans l'est de l'Ile-de-France (3E,48N).
\label{fg:Neige}}
\end{figure}

On a considéré jusque là des points tropicaux, là aussi pour une bonne raison : les premiers tests effectués sur un point en France m'ont donné du fil à retordre avant de comprendre que le problème venait de la neige, puis du fil à retordre ensuite pour y apporter une solution.
En gros, en présence de couvert neigeux, le sol est isolé de l'atmosphère, et les flux radiatifs sont fortement modifiés.
Je n'ai pas trouvé de solution élégante pour prendre en compte la neige dans ces calculs.
En revanche, j'ai inclus une pondération des flux par $1-f_{sn}$ (fraction de neige disponible dans les archive CMIP) et en ajoutant un terme de rappel très fort vers la température de surface lue dans l'archive, de la forme $f_{sn}*\lambda ( \tslu - \ts )$.
Avec cette astuce, on est capable d'obtenir des résultats aussi bons sur les points avec neige que sur les autres (\fig{Neige}).
Evidemment, on "triche" quand il y a de la neige.
Mais ca permet de vérifier que le calcul hors période de neige est valable là aussi.

\subsection{Différence de température entre le bitume et une surface naturelle}


\begin{figure}
% series.sh suivi de inertie_post.sh
\centerline{\includegraphics[angle=-90,width=13cm]{cc.pdf}}
\caption{Comparaison des températures moyennes d'août (noir), température maximum d'août (rouge) et température maximum calculée avec les mêmes variables d'entrée mais en prenant un albédo 0.1 (au lieu de 0.17) et sans prendre en compte l'évaporation (vert). La figure du bas montre l'évolution de la différence de température maximum sur bitume et sur surface naturelle (noir) et du coefficient d'aridité $\beta$. On prend ici le scenario SSP585, le plus pessimiste. Localisation dans l'est de l'Ile-de-France (3E,48N).
\label{fg:Neige}}
\end{figure}

Pour illustration de l'utilisation de ce type de modèle pour des calculs de type "services climatique",
on regarde l'effet sur le calcul des températures maximum d'août 
du fait de considérer une surface de bitume (avec un albédo de 0.1 au lieu de 0.17 et sans évaporation) plutôt que la surface naturelle (celle de la maille du modèle de climat).
La surface de bitume est typiquement une vingtaine de degrés plus chaude que la surface de la maille du modèle.
On voit que cet excès de température sur la surface bituminée tend à décroître légèrement au fur et à mesure que le climat se réchauffe, ce qu'on peut relier au fait que le coefficient d'aridité décroît avec le réchauffement. On voit clairement que dans ce scenario SSP585 (le plus pessimiste), on finit par être en canicule/sécheresse tous les étés, ce qui diminue en partie l'impact de la différence entre la surface naturelle et le bitume.
Ce calcul est sans doute relativement réaliste si on peut supposer que la température de l'atmosphère et le coefficient d'échange turbulent ne sont pas affectés par la surface de bitume, ce qui sera typiquement vrai pour une surface suffisamment petite (une route par exemple).

%-------------------------------------------------------------------------------
\section{Intégration par méthode de Monte Carlo}
%-------------------------------------------------------------------------------

\section{Monte-Carlo de la version discrétisée}

A travailler à Roffiac ?


\section{Sègment $\delta$}

\section{Marche sur sègments}


%-------------------------------------------------------------------------------
\section{Conclusion}
%-------------------------------------------------------------------------------

On a un cadre de validation des modèles d'installation thermique.

On pourrait coder les Hyp1/2/3 à la fois en explicite et implicite.

Je n'ai montré que des figures pour le début des simulations historiques (10 ans) mais j'ai déjà les fichiers météo sur toute la période.


\end{document}
